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SUMMARY 
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results  of  the  use  of  this  code,  also  named  FOURFIT,  and  a  companion  code, 
FOURPLT,  which  permits  the  results  to  be  plotted.  Fits  to  record  traces 
from  two  separate  simulation  events  are  compared  to  previously  published 
results  which  were  determined  using  the  graphical  version  of  the  technique. 
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SECTION  1 


INTRODUCTION 

The  development  and  application  of  a  Fourier  domain  technique  for 
estimating  the  equivalent  nuclear  yield  (W)  and  peak  overpressure  (P  ) 
of  airblast  simulation  records  are  presented  in  References  1  and  2.  The 
estimates  are  made  graphically  by  comparing  the  Fourier  amplitude  spectra 
of  the  data  records  to  the  amplitude  spectra  of  ideal  nuclear  airblast 
curves.  Additionally,  the  methodology  includes  a  means  for  identifying, 
through  low  pass  filtering  of  the  data,  a  "fidelity  frequency"  below  which 
the  data  and  the  ideal  curves  are  in  good  agreement. 

The  technique,  called  FOURFIT,  has  several  advantages.  These 
i ncl ude: 

•  The  Fourier  amplitude  representation  of  the  data 
provides  considerable  insight  into  the  frequency 
content  of  the  data. 

•The  method  provides  consistent  results  for  the  values 
of  W  and  Pso  estimated  for  multiple  records  from 
the  same  event.  Other  methods  give  more  scattered 
estimates  (see  Ref.  3). 

•  The  technique  is  quick  and  easy  to  perform. 

•  FOURFIT  can  be  performed  graphically. 

Despite  the  advantage  of  the  physical  insight  provided  to  an  analyst 
through  graphical  fitting,  an  alternative  "automatic"  fitting  approach  is 
desirable  to  enable  quicker  turn-around  time  for  the  results  of  large 
numbers  of  data  records. 


This  report  documents  such  a  method.  A  computer  program,  FOURFIT, 
has  been  written  which  seeks  to  minimize  the  sum  of  the  squares  of  the 
difference  between  the  data  Fourier  amplitude  and  the  amplitudes  of 
candidate  ideal  nuclear  fits.  The  amplitudes  for  the  candidate  fits  are 
determined  by  an  equation,  parametic  in  P$0  and  W,  which  describes  the 
spectra  for  the  Speicher-Brode  nuclear  overpressure  (Ref.  4).  In 
addition,  the  program  provides  an  estimate  of  the  low  pass  fidelity 
frequency  mentioned  above. 

Section  2  of  this  report  reviews  the  background  of  the  FOURFIT 
technique.  Section  3  discusses  the  structure  and  use  of  the  computer  code 
and  its  algorithm  for  determining  best  fit  equivalent  yield  and  pressure. 
Section  4  presents  some  initial  results  determined  by  the  code.  Section  5 
discusses  some  considerations  into  the  effects  of  high  pass  and  band  pass 
filtering  of  airblast  simulation  data.  Finally,  Section  6  presents 
conclusions  and  recommendations. 


SECTION  2 


THE  FOURFIT  TECHNIQUE 

2.1.  REASONS  FOR  FOURFIT  ANALYSIS 

High  explosive  simulations  of  nuclear  airblast  often  lead  to 
inherent  differences  between  the  simulated  environment  and  the  ideal 
nuclear  environment  being  modeled.  For  example,  high  frequency  spikes  in 
the  early  time  portion  of  High  Explosive  Simulation  Technique  (HEST) 
simulated  airblast  records,  as  seen  in  Figure  1,  are  not  normally  present 
in  actual  nuclear  overpressure  pulses. 

High  frequency  spikes  in  the  HEST  waveform  and  other  simulation 
pressure  history  differences  pose  difficulties  when  fitting  these  records 
with  an  ideal  nuclear  pulse  in  the  time  domain.  Time  domain  fitting 
methods  and,  indeed,  the  acceptance  of  HEST  as  a  useful  simulator,  assume 
that  the  high  frequencies  of  the  HEST  waveform  do  not  drive  the  response 
of  systems  of  interest.  Yet  time  domain  fitting  is  complicated  by  the 
fact  that  high  frequencies  and  low  frequencies  are  superimposed  in  that 
domain.  This  complication  is  especially  difficult  in  the  interpretation 
of  HEST  peak  overpressure  because  the  absolute  peak  overpressure  is 
associated  with  a  high  frequency  spike. 

However,  when  viewed  in  the  frequency  domain,  the  relative 
importance  of  waveform  differences  is  revealed.  The  Fourier  transform 
unfolds  the  various  frequency  contributions  to  the  pressure  history  and 
allows  the  analyst  to  fit  those  spectral  portions  of  a  record  which 
dominate  the  power  in  the  waveform.  Use  of  the  Fourier  amplitude  spectrum 
for  fitting  purposes  thus  provides  for  a  more  accurate  ideal  nuclear  fit 
to  the  simulation  data  than  do  time  domain  techniques. 


2.2.  DEVELOPMENT  OF  THE  FOURFIT  TECHNIQUE 


The  nuclear  airblast  pressure  waveform  satisfies  the  conditions  for 
existence  of  the  Fourier  integral  transform.  That  is, 

•  It  contains  a  finite  number  of  minima  and  maxima. 

•  It  contains  a  finite  number  of  discontinuities. 

•  The  function  is  aperiodic. 

Therefore,  a  measured  airblast  waveform  may  be  Fourier  transformed  using  a 
fast  Fourier  transform  (FFT)  based  upon  the  integral  discrete  Fourier 
transform  (DFT)  as  represented  by  equation  1  below. 

H  ( n/Nat)  =  T  V1  h  (k4t)e'J2,,kn/N  (II 

c  k=0  c 

where  T  =  duration  of  the  signal 

At  =  timestep  between  data  points 

k,n  =  integer  values  and  represent  the  periodicity  of  the  time  and 
the  frequency  functions,  respectively. 

The  DFT  represents  a  limited  duration  signal,  h,  as  one  period  of  an 
infinite  periodic  series  summed  over  N  samples  of  data  and  the  subscript  c 
above  is  used  to  denote  an  approximation  caused  by  this  truncation  of  the 
signal.  The  FFT  computes  a  real  portion  and  an  imaginary  portion  of  the 
Fourier  transform,  H.  These  portions,  in  turn,  may  be  used  to  compute 
Fourier  amplitude  and  phase.  Figure  2  shows  the  Fourier  amplitude 
representation  of  the  pressure  history  for  the  HEST  record  of  Figure  1. 

To  determine  the  nuclear  representation  of  the  simulation  data,  the 
Fourier  spectrum  computed  for  that  data  must  be  compared  to  the  spectra 
computed  for  ideal  nuclear  waveforms.  The  studies  of  References  1  and  2 
and  the  example  presented  in  this  section  were  based  upon  the  "New  Brode" 
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description  of  the  nuclear  waveform  (Ref.  5).  However,  a  later  study 
(Ref.  3)  and  the  FOURFIT  code  discussed  in  following  sections  of  this 
report  are  based  upon  the  more  recent  Speicher-Brode  formulations 
(Ref.  4).  Figure  3  presents  the  "New  Brode"  pressure  histories  for 
several  overpressures  in  scaled  form.  The  scalars  make  the  curves 
applicable  for  all  yields. 

2.2.1.  Estimation  cr  Peak  Pressure  and  Yield 

A  parametric  review  of  Fourier  amplitude  spectra  and  dimensional 
consideration  of  the  DFT  and  the  Brode  equations  revealed  the  proper 
scaling  parameters  for  the  airblast  Fourier  amplitude  spectra.  It  was 
found  that  Fourier  amplitude  scales  by  the  product  of  peak  overpressure 
and  the  cube  root  of  the  yield  (i.e.,  (Pso  x  W1//3)_1).  Furthermore, 
Fourier  frequency  was  found  to  scale  by  the  cube  root  of  the  yield  (i.e., 
W1/3).  Figure  4  illustrates  a  set  of  surface  burst  Brode  Fourier 
amplitude  spectra  in  scaled  form. 

The  normal ization  of  the  ideal  nuclear  spectra  provides  the 
basis  for  the  implementation  of  the  FOURFIT  technique.  The  analyst  first 
notes  that  the  slope  of  each  Brode  spectrum  is  inversely  proportional  to 
the  value  of  peak  overpressure.  That  is,  c.i  Pso  increases  the  spectra 
flatten  as  the  amount  of  power  carried  within  the  higher  frequencies 
increases  relative  to  the  power  in  the  low  frequency  regime.  With  this  in 
mind,  the  analyst  makes  an  initial  estimate  of  the  data  equivalent  peak 
pressure  by  comparing  the  data  spectrum  to  the  Brode  spectra.  The  ideal 
spectrum  which  best  compares  to  the  data  defines  the  initial  estimate  of 

Pso.  The  corresponding  yield  is  determined  as  guided  by  the  scaled 
amplitude  plots.  The  steps  discussed  below  and  illustrated  by  Figure  5 


utilize  those  scaled  plots  as  an  overlay  to  the  data  to  find  the  final 
fit.  The  technique  will  be  illustrated  using  the  HEST  spectrum  shown  in 
Figure  2. 

In  the  first  step,  a  comparison  of  the  data  spectrum  to  the 
Brode  spectra  indicates  a  peak  pressure  of  about  3  megaPascal s  to  be  a 
reasonable  estimate.  Next,  note  that  since  amplitude  scales  by 
(Pso  x  W1/3)-1  and  frequency  scales  by  the  amplitude  scalar 

is  exactly  PSQ  times  the  frequency  scalar.  With  P$0  already 
estimated,  the  problem  is  reduced  to  finding  one  unknown,  namely,  the 
yield.  Graphically,  this  is  performed  by  overlaying  the  spectral  data 
onto  the  ideal  spectra  such  that,  at  equal  locations  on  the  respective 
frequency  axes 

i.e.,  Fb  x  W1/3  =  Fd  (2) 

the  amplitude  axes  are  shifted  by  a  ratio  equal  to  the  estimate  for  Psq 

1-e-.  Pso  x  (Ae/(PS0  x  «1/3))  =  Ad  (3) 

or,  in  the  example,  a  vertical  shift  by  a  factor  of  about  3.  The 
subscripts  B  and  D  refer  to  the  Brode  and  the  data,  respectively.  Since 
the  plots  are  in  the  logarithmic  domain  any  shift  represents  a  constant 
mu 1  tipi ier. 

1  /3 

Similarly,  the  scalar  W  ,  as  a  constant  multiplier,  can  be 
represented  by  a  shift.  This  shift  is  defined  such  that,  since  is 
the  same  multiplier  on  both  axes,  only  a  shift  along  a  line  of  slope  equal 
to  1:1  will  maintain  an  equal  axis  shift.  Furthermore,  the  shift  occurs 
along  a  line  of  -1  slope  since  the  actual  scalars  are  (wl/^j-l  for 
amplitude  and  W^3  for  frequency.  With  the  shift  thus  defined,  it 
simply  is  left  to  perform  this  shift  until  the  data  spectrum  properly 


interpolates  the  Brode  spectra  to  a  location  representing  the  value  for 
peak  overpressure  chosen  for  the  fit  (i.e.,  between  the  2  MPa  spectrum  and 
the  5  MPa  spectrum).  The  resulting  pressure/yield  pair  is  then  computed 
by  "unsealing"  the  overlay.  That  is,  the  comparisons  of  the  respective 
frequency  and  amplitude  axes  for  the  data  and  for  the  normalized  Brodes, 
as  shown  in  the  equations  below,  will  result  in  estimates  of  the  two 
variables  (Pso  and  W). 


Fb  X  w 


Ff  (Wl/3rl 
7TJ7J  ~  (wf  } 


f  _  .  P  (  rf 

CAB/(Pso  X  wl/3)]  S°Vb  X  wl 


The  subscript  f  refers  to  the  fit  to  the  data. 

The  results  for  the  example  are  shown  on  Figure  5  and  the  yield 
i  s  computed  bel  ow. 

3300  Hz  KT1/3  ,,1/3 

200  Hz  =  wf 

or  (6) 

Wf  =  (1.64  KT1/3)3 

=  4.44  KT 

Also,  the  calculation  below  confirms  the  amplitude  scale  using  the  peak 
pressure  (2.78  MPa)  and  the  yield  (4.44  KT)  defined  for  the  fit. 

-1  MP.a,-.s.^  =  (2.78  MPa)  (4.44  KT) 1/3  (7) 

,0218  sec/KTi/J 


A  Brode  overpressure  history  and  associated  impulse  history  and 
Fourier  amplitude  spectrum  are  then  calculated  for  these  values  of  peak 


AT1. 


overpressure  and  yield.  The  three  representations  of  this  signal  can  now 
be  compared  to  those  of  the  data  to  refine  the  fit.  A  need  for  slight 
adjustment  to  the  fit  is  expected  considering  the  graphical  nature  of  the 
technique. 

The  amplitude  spectrum  fit  procedure  is  repeated  to  improve  the 
frequency  and  time  fits.  The  adjusted  estimates  for  the  example  define 
the  pressure  to  be  2.95  MPa  and  the  yield  to  be  5.05  KT.  Figures  6  to  8 
show  the  pressure  and  impulse  histories  and  the  amplitude  spectrum  of  the 
Brode  defined  by  these  final  values  compared  to  the  respective 
representations  of  the  record.  This  is  an  acceptable  fit  and  thus  will  be 
carried  to  the  next  step  of  the  process. 

2.2.2.  Estimation  of  Fidelity  Frequency 

It  may  be  noted  from  Figure  4  that  the  low  frequency  regime  of 
the  nuclear  airblast  carries  significantly  more  power  than  is  carried  by 
the  high  frequencies.  Furthermore,  most  systems  of  interest  in  airblast 
simulations  respond  to  low  frequency  input.  This  section  describes  a 
methodology  which,  through  low  pass  filtering,  quanti tati vely  defines  a 
frequency  cutoff  where  systems  sensitive  to  frequencies  below  that  cutoff 
experience  the  loading  as  defined  by  the  FOURFIT  results  in  the  preceding 
paragraphs.  This  frequency  is  identified  as  the  fidelity  frequency  for 
the  record  since  it  defines  a  limit  above  which  the  fidelity  of  the 
simulation  relative  to  the  prototype  degrades. 

Identification  of  a  fidelity  frequency  makes  use  of  the  fact 
that  most  of  the  power  in  the  airblast  signal  exists  in  the  low 
frequencies.  This  suggests  that  although  low  pass  filtering  of  the  signal 
may  alter  the  exact  shape  of  the  waveform,  it  should  have  limited  effect 


on  the  total  power  delivered.  This  is  illustrated  in  Figures  9  and  10  for 
a  Brode  pulse.  In  Figure  9  it  is  seen  that  as  the  cutoff  frequency  is 
lowered,  the  peak  of  the  filtered  pressure  history  decreases.  However, 
this  peak  reduction  is  accompanied  by  a  broadening  of  that  peak.  This 
suggests  that  the  deliverance  of  power  is  only  delayed,  not  diminished. 
This  hypothesis  is  supported  by  noting  that  the  total  impulse  of  the 
filtered  Brode  (Fig.  10)  is  virtually  unchanged. 

The  determination  of  fidelity  frequency  utilizes  the 
relationships  between  peak  attenuation  and  cutoff  frequency.  It  was  found 
that  the  filtered  peaks  of  the  data  record  show  a  similar  trend  versus  the 
respective  cutoff  frequency.  Fidelity  frequency  is  identified  by  using 
this  similarity  and  with  an  overlay  technique  similar  to  that  used  to 
determine  P$o  and  W  using  the  Fourier  amplitude  spectra. 

A  plot  of  the  data  peak  attenuation  versus  cutoff  frequency  is 

an  overlay  to  the  Brode  peak  attenuation  curves  (Fig.  11).  Note  that  the 

filtered  Brode  peak  is  scaled  by  Psq  and  that  the  Brode  cutoff  frequency 
1  /3 

is  scaled  by  W  .  As  illustrated  in  Figure  11  for  the  record  of  the 
example,  the  data  attenuation  plot  is  shifted  the  proper  amount  on  the 
vertical  axis  and  on  the  horizontal  axis  as  defined  by  Pso  and  by  W, 
respectively.  This  allows  that  the  data  attenuation  curve  approaches  and 
eventually  merges  into  the  attenuation  curve  for  the  equivalent  value  of 
P$0.  The  merger  of  these  curves  defines  the  fidelity  cutoff  frequency 
(180  Hz  for  the  example).  Below  this  frequency  the  data  record  is  similar 
to  the  equivalent  Brode.  Above  this  frequency,  the  data  contains 
diversions  from  the  Brode  which,  in  effect,  do  not  load  systems  sensitive 
only  to  frequencies  below  this  cutoff.  Figure  12  compares  the  filtered 


SECTION  3 


PROGRAM  FOURFIT 

Despite  the  advantages  of  performing  the  graphical  FOURFIT  technique 
described  earlier,  a  quick  "automated"  version  of  the  method  was  desired 
to  provide  rapid  analysis  of  the  typical  airblast  simulation  record. 
Program  FOURFIT  was  written  to  achieve  this  goal.  Appendix  A  provides  a 
1 isting  of  this  code. 

Program  FOURFIT  finds  the  best  fit  nuclear  waveform  for  airblast 
data  based  upon  a  search  to  minimize  the  sum  of  the  squares  of  the 
differences  between  the  data  Fourier  amplitude  spectrum  and  trial  ideal 
spectra.  The  Fourier  transform  of  the  data  is  performed  by  a  call  to  the 
International  Mathematics  and  Statistics  Library  (IMSL)  routine  FFTRC 
(Ref.  6).  IMSL  is  maintained  on  the  Defense  Nuclear  Agency  ( DNA)  CDC 
Cyber  176  computer  on  which  FOURFIT  was  written.  The  ideal  nuclear 
Fourier  amplitudes  for  trial  fits  are  computed  by  an  equation,  parametric 
in  Pso  and  W,  which  describes  the  amplitude  spectra  for  a  set  of  curves 
similar  to  those  presented  in  Figure  4.  The  equation,  however,  describes 
the  more  recent  Speicher-Brode  fit  to  the  nuclear  overpressure  waveform 
(Ref.  4)  represented  in  scaled  form  in  Figure  13.  Scaled  Speicher-Brode 
impulse  histories  and  Fourier  amplitude  spectra  are  presented  in  Figures 
14  and  15,  respectively.  A  detailed  description  of  the  fitting  program 
fol 1 ows. 

Thp  program  is  run  by  reading  an  input  deck  (TAPE2)  to  define  user 
options.  Results  are  listed  in  output  (TAPE6)  and  plotting  information  is 
written  to  TAPE48  to  be  plotted  by  FOURPLT,  the  accompanying  plotting 


program. 


3.1.  INPUT  VARIABLES 


The  main  purpose  of  the  program  FOURFIT  is  to  perform  fits  to 
surface  burst  airblast  simulation  data.  The  program  compares  data  Fourier 
amplitude  spectra  to  estimates  of  the  Fourier  amplitude  spectra  of  the 
ideal  nuclear  represented  by  the  Speicher-Brode  pressure-time  equation. 
However,  the  code  is  also  set  up  to  perform  FFT  analysis  of  the  data 
without  finding  a  fit,  or  to  perform  a  FFT  on  a  specific  Speicher-Brode 
defined  by  a  pressure/yiel d  pair  requested  by  the  analyst.  FOURFIT  also 
computes  the  impulse  history  for  either  of  these  latter  cases.  Finally, 
the  code  is  set  up  to  perform  a  low  pass  Butterworth  filter  of  either  data 
or  a  Speicher-Brode  at  up  to  seven  cutoff  frequencies,  specified  by  the 
user,  per  run. 

The  options  mentioned  above  are  to  be  chosen  by  the  analyst  and  read 
from  an  input  deck,  TAPE2,  assembled  by  that  analyst.  The  contents  of 
TAPE 2  are  summarized  below.  The  code  reads  all  input  lines,  including 
those  not  used  in  analysis,  in  all  runs. 

Line  1:  NEPTS,  IUNITS,  JUNITS  (Format  315) 

The  value  for  NEPTS  is  the  number  of  points  to  be  read  from  the  data 
tape.  The  other  variables  in  this  line  account  for  data  input  units.  For 
IUNITS  greater  than  zero,  the  code  assumes  that  the  pressure  values  are 
read  in  pounds  per  square  inch  and  converts  the  data  to  megaPascal s,  the 
internal  units  of  the  code.  For  IUNITS  less  than  zero,  the  data  is 
assumed  to  be  read  in  the  program  units.  Furthermore,  the  program  works 
in  units  of  seconds  for  time  and  Hertz  for  frequency.  JUNITS  less  than 
zero  indicates  an  input  time  step  consistent  with  this  fact.  For  JUNITS 
greater  than  zero,  the  program  assumes  input  in  milliseconds  and  performs 


the  proper  conversion.  The  value  of  these  variables  do  not  affect  the 
calculation  of  a  Spei cher-Brode  function. 

Line  2:  PSOI,  WI  (Format  2F5.2) 

The  meaning  of  these  terms  in  line  2  differs  depending  on  the 
program  option  (explained  in  line  3)  chosen.  For  analysis  of  a  data  trace 
for  its  impulse  and  FFT  or  for  filtering  that  data,  without  performing  a 
fit,  PSOI  and  WI  are  not  used.  However,  for  the  case  of  fitting  the  data 
with  a  Spei cher-Brode,  these  variables  provide  the  "seeds"  for  defining 
the  candidate  fits.  PSOI  is  the  seed  for  peak  overpressure  in  MPa  and  WI 
is  the  seed  for  nuclear  yield  in  kilotons.  A  wide  range  about  these  seeds 
(plus  and  minus  a  decade  for  each)  is  tested  and  so  they  need  not  be 
exceptionally  close  to  the  final  values.  However,  a  good  set  of  seeds  may 
nominally  be  considered  to  be  the  event  design  pressure  and  yield.  For 
cases  in  which  only  the  Spei cher-Brode  will  be  analyzed,  PSOI  and  WI  are 
the  peak  overpressure  and  yield  values,  respectively,  of  the  ideal 
waveform  to  be  calculated. 

Line  3:  IOPT,  IFILT  (Format  215) 

IOPT  defines  the  program  option  to  be  run.  IOPT  equals  1  for 
performing  the  FOURFIT  automated  fitting  routine.  To  simply  integrate  and 
FFT  the  data,  IOPT  equals  2.  For  IOPT  equals  3,  the  code  analyzes  only 
the  Speicher-Brode  specified  by  PSOI  and  WI  in  line  2.  The  value  for 
IFILT  determines  whether  or  not  the  pressure  history  is  to  be  filtered. 

No  filtering  is  done  for  IOPT  equal  to  2  or  3  if  the  value  of  IFILT  is 
less  than  zero.  A  value  of  IFILT  greater  than  zero  for  these  same  IOPT 
values  performs  a  low  pass  Butterworth  filter  on  the  pressure  history  at 
the  cutoff  frequencies  defined  by  FLO  (line  4).  For  IFILT  greater  than 
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zero,  neither  impulse  nor  FFT  calculations  are  performed.  For  IOPT  equal 
to  1,  IFILT  may  be  any  value. 

Line  4:  FLO(I),  1=1,7  (Format  7F10.0) 

FLO(I)  are  the  low  pass  cutoff  frequencies  in  Hertz  used  by  the 
Butterworth  filter.  Up  to  seven  filter  levels,  in  any  order,  are  allowed 
for  options  2  and  3.  For  the  number  of  filters,  N,  less  than  seven, 

FL0(N  +  1)  must  be  set  to  0.  Although  option  1  does  not  utilize  the  IFILT 
value  (line  3),  it  nevertheless  performs  filtering  in  order  to  determine 
fidelity  frequency.  Therefore,  IOPT  =  1  requires  that  TAPE2  contain 
several  filter  levels  to  be  defined  on  this  line.  Furthermore,  the 
algorithm  requires  that  these  filters  be  in  descending  order 
(e.g.,  FL0(1)  =  5000.,  FL0(2)  =  2000.,  FL0(3)  =  1000.,  etc.).  The  values 
for  FLO  in  this  series  may  be  chosen  by  the  analyst.  However,  bounding 
values  of  FL0(1)  =  5000.  and  FL0(7)  =  50.  with  a  reasonable  spread  of 
values  for  FLO(I),  I  =  2,6  between  these,  have  proven  to  be  adequate. 
Additionally,  each  FLO(I)  must  be  less  than  or  equal  to  the  Nyquist 
frequency  for  the  digital  filter  to  remain  stable.  (Note  that  if  either 
FL0(1)  orFL0(7),  i.e.,  the  limiting  cases,  is  determined  to  be  the 
fidelity  frequency,  the  values  should  be  altered  accordingly  and  the 
program  resubmitted.) 

These  four  lines  complete  the  input  deck  needed  to  submit  a  FOURFIT 
run.  An  example  input  deck  is  listed  below.  The  program,  will 
subsequently  compute  a  6000  point  FFT  on  data  read  in  the  units  of  psi  and 
seconds  (line  1).  The  initial  estimated  pressure  and  yield  pair  are 
20.  MPa  and  2.  KT,  respectively  (line  2).  The  third  line  indicates  that  a 
fit  will  be  performed.  The  second  value  in  this  line  represents  the 


filter  switch  and,  since  this  run  asks  for  a  fit,  it  will  not  be  used.  The 
final  line  of  input  lists  the  seven  low  pass  filter  levels,  in  Hz,  to  be 
tested  for  locating  the  low  pass  fidelity  frequency. 

Column 

5  10  15  20  30  40  50  60  70 

6000  1  -1 

20.  2. 

1  -1 

5000.  2000.  1000.  500.  200.  100.  50. 

3.2.  PROGRAM  STRUCTURE 
3.2.1.  Data  Calculations 

IOPT  equal  to  2  is  the  simplest  option  to  perform  and,  hence,  the 
option  taking  the  most  direct  calculation  route.  This  option  simply  requires 
a  subroutine  to  read  the  data,  a  subroutine  to  integrate  that  data  for 
impulse  and  a  subroutine  to  calculate  the  Fourier  transform  of  the  record. 

For  IFILT  greater  than  1  (filters  to  be  executed)  the  impulse  and  FFT  are  not 
calculated.  Instead,  a  filter  subroutine  is  called  and  the  record  is  low 
pass  filtered  at  specified  cutoff  frequencies. 

The  program  FOURFIT,  as  listed  in  Appendix  A,  calls  the  subroutine 
EBREAD  to  read  the  data  pressure  histories.  EBREAD  is  set  up  to  read  a  card 
image  format  (EBCDIC)  tape  of  the  data  and  its  header.  The  format  of  EBREAD 
is  the  format  of  several  tapes  analyzed  by  the  author  which  were  provided  by 
the  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES).  The  format  that 
those  tapes  employed  was  pressure  data  written  as  five  data  values  per  card 
(5E16.8).  The  set  of  cards  for  a  given  trace  is  preceded  by  a  header  record 
containing  shot  and  data  information  (i.e.,  shot  title,  gage  title,  time 
step,  total  number  of  points)  in  the  format  of  3(2A10),  E15.8,  15.  These 
tapes  have  been  written  in  psi  for  pressure  and  the  program  converts  the 
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values  to  MPa.  The  time  step  is  written  in  seconds.  Any  tape  of 
different  format  must  be  accompanied  by  a  substitution  for  EBREAD  to  read 
that  data.  However,  this  new  subroutine  must  retain  the  structure  of 
EBREAD  if  the  program  is  to  perform  properly.  This  includes  proper  units 
(pressure  in  MPa,  time  in  seconds),  proper  ordering  of  calls  to  impulse, 
FFT  and  filtering  subroutines  and  identical  writes  to  TAPE48. 

After  the  data  has  been  read,  IOPT  equal  to  2  causes  EBREAD  to 
take  one  of  two  paths,  depending  on  the  value  of  IFILT.  If  filtering  is 
not  to  be  done,  the  subroutine  causes  the  impulse  history  and  the  Fourier 
amplitude  spectrum  to  be  calculated.  Subroutine  IMPULSE  integrates  the 
data  by  Simpson's  approximation.  Subroutine  FFTRC  computes  the  fast 
Fourier  transform  after  the  algorithm  of  Singleton  (Ref.  7).  On  the  other 
hand,  if  filter  histories  are  desired,  subroutine  FILTER  filters  the  data 
using  the  recursive  equations  derived  for  a  two  pole  low  pass  Butterworth 
filter  as  found  in,  for  example.  Reference  8. 

3.2.2.  Speicher-Brode  Calculations 

IOPT  equal  to  3  causes  calculation  of  a  Speicher-Brode  pressure 
history  and  either  its  impulse  and  FFT  or  its  specified  low  pass  filtered 
pressure  histories.  Structure  of  subroutine  SPBRODE  is  similar  to  that  of 
EBREAD  except  that  the  data  reads  are  substituted  for  by  the 
Speicher-Brode  equations.  Also,  SPBRODE  requires  a  target  range  to 
perform  its  calculations.  Therefore,  before  entering  into  SPBRODE,  the 
program  utilizes  subroutines  RANGE  and  PPEAK  to  iterate  on  the  distance 
from  surface  burst  ground  zero  for  the  given  PSOI  and  WI  pair  specified. 
Impulse  and  FFT  or  filter  histories  are  calculated  as  described  for  IOPT 
equal  to  2. 
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3.2.3.  Fit  to  Data 

A  value  of  IOPT  equal  to  1  causes  th.e  code  to  find  a  best  fit 
nuclear  waveform  for  the  data.  That  data  is  read,  integrated  and  Fourier 
transformed.  A  least  squares  algorithm  finds  a  best  fit  to  the  data 
Fourier  amplitude  spectrum  based  on  an  equation,  parametric  in  peak 
pressure  and  yield,  which  describes  the  scaled  Speicher-Brode  Fourier 
amplitude  spectra.  The  actual  best  fit  Speicher-Brode  waveform  is 
calculated,  integrated  and  Fourier  transformed  for  comparison  to  the  data 
in  the  form  of  plots.  The  data  pressure  history  and  the  equivalent 
Speicher-Brode  are  low  pass  filtered  at  several  levels  of  frequency 
cutoff.  The  peaks  of  these  filtered  histories  are  compared  in  the 
determination  of  fidelity  frequency. 

The  best  fit  search  is  performed  within  the  subroutine  FIT. 

This  subroutine  is  modeled  after  a  similar  subroutine  discussed  in 
Reference  9.  The  search  begins  with  a  set  of  five  peak  pressure  values 
and  five  yield  values.  These  values  are  equivalent  to  the  product  of  the 
coefficients  .1,  .4,  1.,  4.  and  10.  times  PSOI  and  WI.  The  final  results 
for  equivalent  peak  overpressure  and  yield  are  found  within  these  values. 
(If  the  final  result  for  either  pressure  or  yield  is  either  .1  or  10. 
times  PSOI  or  WI,  respectively,  i.e.,  the  limiting  cases,  the  analyst  is 
advised  to  alter  the  seeds  accordingly  and/or  to  check  the  quality  of  the 
data. ) 

Subroutine  FIT  takes  a  value  of  peak  overpressure,  PP,  and  pairs 
that  value  with  each  of  the  five  yield  values,  W.  The  comparison  between 
each  candidate  Speicher-Brode  and  the  data  occurs  after,  for  each  value  of 
frequency  for  the  data  spectrum,  a  value  of  ideal  nuclear  Fourier 
amplitude  is  computed  using  the  parametric  equation  below. 


A2  =  .  01474PP 


W  ' 


A3  =  .0011PP 


(PP)'-23YFs  \-2.15 


A4  =  . 001 32( F  ) 


-.547 


A5  =  . 01034PP- *113(Fs) -1* 


A6  =  . 00001 1PP 


.77/  Fs 


. 3/  Fs 


A7  =  . 0000666PP 


ASCL  *  A1  -  A2  +  A3  +  A4  ♦  A5  -  A6  +  A7  (8h) 

A$b  =  ASCL  x  PP  x  W1/3  (81) 

where  F$  =  yield  scaled  frequency  (F  x  W^3) 

Fso  =  yie1d  scaled  fundamental  frequency  (FQ  x  W  '  ) 

FQ  =  1/positive  phase  duration 
ASCL  =  scaled  Spei cher-Brode  Fourier  amplitude 
(A3B/(Pso  x  W1/3„ 

A^g  =  estimated  Speicher-Brode  Fourier  amplitude 
This  equation  represents  a  fit  to  the  normalized  surface  burst  Speicher- 
Brode  Fourier  amplitude  spectra  shown  in  Figure  15.  (Use  of  this  equation 
is  facilitated  when,  for  each  candidate  fit,  successive  calls  to  RANGE  and 
PPEAK  calculate  the  positive  phase  duration.) 


For  a  given  data  frequency,  equation  3  is  used  to  calculate  a 
Speicher-Brode  amplitude  for  the  trial  peak  pressure  and  yield  pair.  In 
order  to  assimilate  the  graphical  methodology,  which  is  carried  out  in 
log-log  form,  the  common  logarithm  of  the  amplitude  estimate  is  computed. 
The  difference  between  this  last  value  and  the  common  logarithm  of  the 
data  amplitude  is  computed  and  then  is  squared  to  accomodate  algebraic 
sign.  To  further  assimilate  the  graphical  method,  this  difference  is 
divided  by  the  particular  data  frequency.  This  last  step  considers  that 
the  FFT  is  computed  using  a  constant  frequency  step.  Therefore,  the  point 
density  increases  by  an  order  of  magnitude  with  each  decade  increase  in 
frequency.  Division  by  the  frequency  compensates  for  the  weighting  that 
resul ts. 

The  value  just  computed  is  added  in  a  summation  of  squared 
differences  for  the  PP/W  pair  for  the  present  trial.  This  summation 
process  is  repeated  starting  at  the  data  fundamental  frequency  and 
continuing  to  some  high  frequency  (a  value  which  is  a  function  of  the  data 
record)  which  sufficiently  includes  the  signficant  portion  of  the  data. 
(This  final  frequency  must  be  chosen  so  as  to  include  the  peak  of  the 
data,  but  must  be  limited  to  avoid  extensive  calculation  within  the  high 
point  density  "record  noise"  frequency  regime.)  The  results  of  this 
summation  provide  a  value  DELTAW  for  the  present  pressure/yield  pair. 

A  DELTAW  value  is  computed  for  each  of  the  other  candidate 
yields  until,  for  the  given  PP(J),  five  DELTAW(I)  values  are  available  for 
comparison.  The  minimum  of  these  five  values  is  then  located  and  the  next 
iteration  occurs  with  a  new  set  of  yields  replacing  the  old  set  with  the 
yield  which  gave  the  minimum  DELTAW(I)  being  the  central  value.  (For 
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example,  for  W{I),  I  =  1,5,  DELTAW(2)  may  have  been  a  minimum.  For  the 
next  iteration,  then,  W ( 2 )  becomes  the  new  W(3);  W(l)  remains  unchanged 
and  the  past  W( 3)  becomes  the  new  W(5).  The  new  W(2)  and  W(4)  values  are 
intermediate  to  the  new  W(l)  and  W ( 3 )  and  the  new  W(3)  and  W(5), 
respectively.)  This  procedure  is  repeated  until,  for  the  given  PP(J),  the 
field  of  W { I )  is  narrowed  so  that  the  extremes  are  within  1  percent  of 
each  other  (i.e.,  2  x  (W(5)  -  W(1))/(W(5)  +  W(l))  £  .01).  When  this 
tolerance  is  met,  the  minimum  DELTAW(I)  of  the  final  group  is  set  to  be 
DELTAP(J)  for  the  given  value  of  PP(J). 

The  iteration  next  proceeds  to  the  second  value  for  PP(J) 
coupled  with  each  of  the  original  five  values  for  yield.  Eventually,  four 
more  DELTAP(J)  values  are  computed  and  the  minimum  of  the  five  DELTAP(J) 
values  is  determined.  In  this  manner*  the  field  of  PP(J),  J  =  1,5  is 
narrowed  to  five  new  values  and  the  entire  process  continues  until  the 
spread  of  PP(0)  values  are  limited  to  within  a  tolerance  as  was  specified 
for  the  W ( I )  above.  When  this  criterion  is  met,  the  best  fit 
Speicher-Brode  is  considered  to  be  found  as  the  pressure/yield  pair  for 
which  the  final  minimum  DELTAP(J)  was  found.  (Recall,  though,  that  if  the 
final  pressure  or  yield  is  either  .1  times  or  10.  times  the  respective 
seed,  the  validity  of  the  answer  must  be  checked.) 

Next,  for  this  case  where  I0PT  =  1 ,  a  Speicher-Brode  pressure 
history  is  computed  for  the  pressure/yield  pair  defined  by  FIT.  This 
history  is  integrated  and  Fourier  transformed  so  that  final  plots, 
provided  by  program  FOURPLT,  represent  comparisons  between  the  data  and 
the  best  fit  Speicher-Brode  pressure  history,  impulse  history  and  Fourier 
amplitude  spectrum.  Additionally,  this  ideal  waveform  is  low  pass 


filtered.  The  peaks  of  the  filtered  pressure  histories  are  determined  and 
compared  to  the  peaks,  at  their  respective  cutoff  frequency,  of  the 
filtered  data  which  have  also  been  determined.  Beginning  with  the  highest 
frequency  filter  and  proceeding  in  descending  order,  the  level  at  which 
the  filtered  data  peak  is  within  10  percent  of  the  filtered  Speicher-Brode 
peak  is  chosen  as  the  low  pass  fidelity  frequency.  If  a  fidelity 
frequency  is  not  found,  a  value  of  -999.  is  assigned  to  this  variable.  If 
this  value  appears  in  the  output,  further  investigation  is  warranted. 

The  analysis  discussed  in  the  previous  paragraphs  creates  a 
file,  TAPE48,  which  contains  information  to  be  plotted.  Program  FOURPLT, 
discussed  below,  utilizes  this  file  to  produce  plotted  output  of  results 
for  the  three  program  options  (IOPT).  In  addition,  the  results  for  IOPT 
equal  1,  fitting  of  data,  will  be  written  to  the  output  file  for  the 
analyst's  record.  An  example  of  this  output  is  shown  in  Figure  16.  (For 
IOPT  equal  to  3,  just  Speicher-Brode  calculations  to  be  done,  similar 
output  is  provided.)  A  flow  chart  of  FOURFIT  is  presented  in  Appendix  B. 
Appendix  C  provides  a  flow  chart  of  subroutine  FIT. 

3.3.  PROGRAM  FOURPLT 

Program  FOURPLT  is  to  be  used  in  conjuction  with  program  FOURFIT. 
FOURPLT  exists  to  attach  the  TAPE48  made  by  FOURFIT,  read  that  file,  known 
as  TAPE9  within  FOURPLT,  and  plot  the  contents.  The  program  uses  the 
DISSPLA  plotting  capabilities  maintained  for  the  DNA  CDC  Cyber  176 


computer.  A  separate  plotting  routine  allows  for  analysis  of  greater 
sizes  of  data  arrays  and  provides  for  quicker  turn  around  by  dividing  the 
core  requirement  of  the  combined  job.  FOURPLT  requires  no  input  other 
than  the  TAPE48  file  to  run  successfully.  A  listing  of  program  FOURPLT  is 
provided  as  Appendix  D. 
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3.4.  PROGRAMMING  NOTES 

Before  continuing  with  presentation  of  the  results  from  running 
FOURFIT  and  FOURPLT,  a  point  should  be  noted  which  will  assist  the 
operator  in  successfully  running  the  codes.  The  Speicher-Brode  Fourier 
amplitude  equations  describe  the  total  positive  phase  duration  of  the 
respective  pressure  histories.  Therefore,  the  data  to  be  analyzed  should 
similarly  be  carried  as  nearly  as  possible  through  full  positive  phase. 

Since  the  amount  of  data  that  can  be  analyzed  is  set  within  several 
arrays,  these  arrays  must  be  large  enough  to  contain  all  of  the  data. 

This  includes  pressure  (PRESS)  and  time  (TTIM)  to  be  dimensioned  at  least 
as  large  as  NEPTS;  impulse  (PIMP)  and  impulse  times  (TIMP)  must  be 
dimensioned  at  least  (NEPTS/2)-l;  amplitude  (AMP)  and  frequency  (FRQ)  must 
be  at  least  as  large  as  (NEPTS/2)+l.  In  addition,  the  data  FFT  working 
arrays  (IWKE,  WKE)  must  be  of  sufficient  size.  (Reference  6  contains  an 
algorithm  for  computing  the  necessary  size  of  these  two  working  arrays  by 
factoring  NEPTS.)  The  Speicher-Brode  calculations  use  these  identical 
arrays  and  the  number  of  points  assigned  to  these  calculations  is  NBPTS 
equal  to  2048.  Each  array  must  be  at  least  large  enough  to  accomodate 
thi s  val ue. 
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SECTION  4 


PROGRAM  RESULTS 

4.1.  SAMPLE  OUTPUT 

Figures  17  to  28  illustrate  the  possible  output  from  FOURFIT  and 
FOURPLT  for  the  various  options.  Figures  17  to  19  result  from  requesting 
IOPT  =  2  with  IFILT  =  -1.  The  program  reads,  integrates  and  Fourier 
transforms  a  data  record,  in  this  case  record  AB-5  from  event  0.35  KBAR 
DISC  HEST.  In  Figure  20,  the  same  data  record  is  read,  IOPT  =  2,  but  is 
low  pass  filtered,  IFILT  =  1,  at  a  cutoff  frequency,  FLO,  equal  to 
1000.  Hz.  In  each  case,  the  identifying  header  is  presented  at  the  top  of 
the  plot. 

In  Figures  21  to  23,  the  results  of  computing,  integrating  and 
Fourier  transforming  a  Spei cher-Brode  pressure  waveform  (PSOI  =  39.60  MPa, 
HI  =  0.87  KT,  IOPT  =  3,  IFILT  =  -1)  are  presented.  Figure  24  presents 
this  same  waveform  {IOPT  =  3)  low  pass  filtered  (IFILT  =  1)  at  1000.  Hz 
(FLO)  frequency  cutoff.  Information  to  the  right  of  the  plot  identifies 
the  ideal  waveform  plotted.  Each  of  the  types  of  submittals  discussed  in 
this  and  the  preceding  paragraph  require  on  the  order  of  0.5  CP  seconds  of 
execution  time  on  the  Cyber  176  computer. 

Finally,  Figures  25  to  27  show  an  example  of  the  results  of  a  run  to 
fit  the  data  record  (IOPT  =  1)  presented  in  Figure  17.  Both  the  data  and 
the  computed  best  fit  Speicher-Brode  are  plotted  for  comparison  of 
pressure  histories,  impulse  histories  and  Fourier  amplitude  spectra. 

Figure  28  compares  the  data  and  its  equivalent  fit  both  low  pass  filtered 
at  the  low  pass  fidelity  frequency  identified  by  the  fit.  The  job 
requiring  a  fit  to  this  record  required  on  the  order  of  80  CP  seconds 


execution  time  on  the  Cyber  176  computer.  (The  required  computer  time 
will  vary  as  the  number  of  data  points  and,  hence,  frequency  comparisons, 
varies.  ) 

4.2.  RESULTS  OF  FITTING  ROUTINE 

Several  records  from  the  test  0.35  KBAR  DISC  HEST  were  fit  using 
program  FOURFIT.  These  records  were  chosen  because  they  were  data 
histories  which  were  fairly  representative  of  a  nuclear  pressure 
waveform.  Furthermore,  these  records  were  previously  analyzed  using  the 
graphical  FOURFIT  technique  (Ref.  3)  and  thus  provided  a  basis  for 
comparison  between  the  program  results  and  accepted  fits. 

The  fits  calculated  for  several  records  from  the  DISC  HEST  are  shown 
in  Figures  29  to  49.  (These  results  are  in  addition  to  those  for  record 
AB-5  shown  earlier.)  These  figures  represent  the  output  from  running 
FOURFIT  and  FOURPLT  and  includes  both  the  data  and  the  fit  compared  by 
pressure  history,  impulse  history  and  Fourier  amplitude  spectrum  for  each 
record.  Table  1  summarizes  these  fits  and  compares  the  results  of  the 
program  ("automated")  to  fits  found  graphically  ( "graphical " ) .  The 
authors  feel  that  Figures  29  to  49  show  favorable  comparisons.  (Though 
some  of  the  latest  automated  results  do  not  agree  closely  with  previous 
graphical  resul ts- -i .e. ,  yield  values  for  AB-4,  for  AB-9  and  for 
AB-10--the  plot  comparisons  for  these  records  are  very  acceptable.) 

The  program  was  next  applied  in  the  analysis  of  a  second  event, 

0.35  KBAR  HEST.  FOURFIT  managed  to  determine  acceptable  fits  to  several 
of  the  records  (e.g.,  the  fit  for  record  51  shown  in  Figures  50  to  52). 
However,  for  other  records  of  this  event,  the  data  record  and  its  time 
domain  fits  diverge  after  several  milliseconds.  This  divergence  may  be 
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traced  to  trends  in  the  pressure  data  which  are  atypical  of  the  ideal 
nuclear  pulse.  These  trends,  as  noted  for  one  of  the  records,  are 
discussed  below. 

Figure  53  presents  the  pressure  history  of  record  417  from  0.35  KBAR 
HEST.  Several  variations  between  the  ideal  waveform  and  record  417  are 
readily  noticeable.  First,  in  contrast  to  the  DISC  HEST  records,  record 
417  contains  one  relatively  high  magnitude  (46  MPa),  but  very  narrow  spike 
suggesting  a  low  amount  of  power  carried  in  the  peak  of  the  record.  This 
spike  is  followed  by  an  extended  vibratory  component  at  about  11  MPa  and 
another  at  about  7  MPa  suggesting  that  a  lower,  but  still  fairly  high 
frequency  regime  may  be  sustaining  too  much  power.  The  remainder  of  the 
waveform  shows,  rather  than  a  purely  exponential  decay,  a  decay  that,  at 
times,  shows  a  nearly  linear  trend,  upon  which  is  superimposed  a  low 
frequency  oscillatory  component.  This  would  foretell  a  rise  in  the 
amplitude  spectrum  in  the  low  frequency  end. 

The  Fourier  amplitude  spectrum  for  record  417  was  computed  and  is 
presented  in  Figure  54.  This  figure  fulfills  the  expectations  resulting 
from  review  of  the  pressure  history.  The  spectrum  falls  off  rapidly 
between  the  fundamental  frequency  and  about  150  Hz.  At  this  point,  the 
slope  changes  to  a  lesser  decay  of  power  toward  the  intermediate  to  high 
frequencies.  As  the  Nyquist  frequency  is  approached,  the  spectrum  falls 
off  more  rapidly.  These  observations  seem  to  correspond  to  the  low 
frequency  oscillation,  the  early  time/low  magnitude  oscillations  and  the 
narrowness  of  the  peak,  respectively. 

The  factors  listed  above  indicate  that  the  recorded  trace  carries  a 
low  qualitative  fidelity  in  comparison  to  the  Speicher-Brode.  These  data 


trends  suggest  that  the  fits  to  such  records  may  show  obvious  variances  in 
comparisons  to  those  records.  Figure  55  presents  the  Fourier  amplitude 
fit  determined  by  a  FOURFIT  run  on  record  417.  It  is  seen  that  the  data 
spectrum  diverges  from  the  fit  at  times,  with  the  non-nuclear  trends  of 
the  data  becoming  obvious.  The  pressure  history  comparison.  Figure  56, 
also  illustrates  regions  of  divergence  between  the  data  and  the  nuclear 
waveform.  The  fit  impulse  history  (Fig.  57)  is  seen  to  be  a  mismatch  to 
the  data  beyond  about  10  msec.  The  difficulties  encountered  in  analysis 
of  poor  fidelity  records,  such  as  417,  indicates  a  need  for  more  study 
into  the  approach  for  analysis  of  such  records.  For  example,  different 
frequency  regimes  of  such  data  may  be  subject  to  varying  weighting 
functions  to  perform  the  fit. 

Comparisons  between  other  records  from  0.35  KBAR  HEST  and  their 
respective  fits  are  shown  in  Figures  58  through  72.  These  fits  are 
summarized  in  Table  2. 


Table  2.  Comparison  of  FOURFIT  results  for  0.35  KBAR 
HEST:  graphical  (ref.  3)  versus  automated. 


SECTION  5 


FILTER  STUDY 

The  previous  sections  discuss  the  development  of  the  FOURFIT 
technique  and  a  computer  code  written  for  the  purpose  of  performing  that 
technique  numerically.  Consideration  of  fidelity  through  low  pass 
filtering  suggested  that  more  information  could  be  culled  from  data 
records  through  more  extended  analysis.  Specifically,  it  was  hypothesized 
that  high  pass  filtering  and  possibly  band  pass  filtering  of  the  data  and 
of  the  nuclear  waveforms  could  prove  insightful.  For  example,  different 
frequency  regimes  of  the  simulated  waveform  may  represent  a  different 
equivalent  peak  pressure  and/or  yield  for  systems  with  sensitivity  in 
those  frequency  regimes. 

An  extended  version  of  FOURFIT  was  written  to  include  high  pass  and 
band  pass  filters.  As  in  the  low  pass  filter,  these  filter  types  were  two 
pole  recursive  digital  Butterworth  filters.  Extensive  analysis  into  the 
effects  of  these  filters  on  the  Speicher-Brode  waveform  was  undertaken.  A 
band  pass  filtered  ideal  waveform  of  Figure  21  is  shown  in  Figure  73.  It 
can  be  seen  that  the  band  pass  filter  drastically  alters  the  form  of  the 
signal.  This  most  likely  is  due  not  only  to  a  removal  of  low  frequency 
power,  but  also  to  phase  shifting. 

Although  low  pass  filtering  left  the  final  airblast  impulse 
virtually  unaffected,  it  is  obvious  from  Figure  73  that  the  impulse  and, 
hence,  the  power  of  the  original  signal  are  reduced.  Several  attempts 
were  made  to  correlate  this  reduction  in  power  to  various  factors. 
Comparisons  using  varying  high  pass  cutoff  and  low  pass  cutoff 
combinations  with  various  P$0/W  pairs  failed  to  produce  any  promising 
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results.  These  comparisons  included  studies  of  the  filtered  peaks  (both 
positive  phase  and  negative  phase  peaks)  and  of  impulse  (both  positive 
phase  and  total  impulses). 

High  pass  filter  studies  were  more  promising  than  those  regarding 
the  band  pass  filter.  Figure  74  presents  a  high  pass  filtered  waveform 
from  Figure  21.  Although  the  waveform  is  altered  similar  to  the  effect  of 
the  band  pass  filter,  a  correlation  between  the  filtered  trace  and  the 
filter  was  developed.  Before  discussing  this  effect,  it  must  first  be 
noted  that  recursive  digital  filters  are  a  function  of  the  data  time  step 
as  well  as  of  the  cutoff  frequency  and  the  number  of  filter  poles.  It  was 
found  that  the  time  step  dependence  is  an  important  factor  for  a  high  pass 
filter  of  a  waveform  such  as  the  nuclear  airblast  pulse  (i.e.,  sudden  rise 
to  peak).  With  this  in  mind,  several  Speicher-Brode  waveforms  of  various 
peak  pressure  and  yield  combinations  were  calculated  with  the  same  time 
step  for  each  and  were  then  filtered  at  various  high  pass  cutoff  levels. 
Figure  75  shows  the  effect  of  cutcff  frequency  on  the  scaled  peak  of  the 
filtered  ideal  waveform.  The  curve  applies  for  all  yields  and  represents 
the  high  pass  filter  peak  attenuation  curve  for  data  with  a  sampling  rate 
of  100  kHz  (the  sampling  rate  of  the  0.35  KBAR  DISC  HEST).  The 
application  of  the  attenuation  curve  to  the  data  analysis  is  described 
bel  ow. 

Given  the  time  step  of  the  data  to  be  analyzed,  a  table  of  values 
for  Speicher-Brode  peak  attenuation  as  a  function  of  those  cutoff 
frequencies  listed  in  the  input  deck  (TAPE5)  must  be  developed.  This 
array  of  information  is  then  added  to  the  code.  Then,  in  the  process  of 
running  FOURFIT,  the  data  must  be  high  pass  filtered  and  the  filtered 


record  peaks  stored.  Next,  the  code  determines  a  best  fit  Speicher-Brode 
and  a  fidelity  frequency.  For  all  systems  with  sensitivities  below  this 
value,  the  equivalent  fit  is  a  valid  loading  function.  Systems  sensitive 
to  frequencies  higher  than  the  fidelity  frequency  will  experience  a 
different  loading  function.  This  different  function  is  determined  by 
reference  to  the  high  pass  attenuation  table  for  the  ratio  of  filter  peak 
to  Pso  for  the  specific  fidelity  frequency  just  determined.  Through 
this  value,  the  high  pass  loading  function  is  determined  by  the  following 
relation: 

(PHPD)FF/(Stored  Rat1'o)FF  =  PSoHP 

or  ^HPD^f/^HPB^soVf  =  PsoHP  (9) 

where  Pj_|pp  =  the  peak  of  the  high  pass  filtered  data 

PHPB  =  Peak  hl9h  Pass  filtered  Speicher-Brode 

psoHP  =  the  high  pass  equivalent  Speicher-Brode  for  the 
data  record. 

The  subscript  FF  refers  to  the  respective  values  at  the  fidelity  frequency. 

No  yield  dependence  was  found  in  this  study.  Therefore,  the  high 
pass  equivalent  waveform  is  assigned  a  yield  identical  to  that  of  the  low 
pass  equivalent  waveform.  Figure  76  presents  a  plot  of  record  AB-5  from 
0.35  KBAR  DISC  HEST  with  its  FOURFIT  comparison.  The  high  pass  equivalent 
peak  overpressure  is  identified  on  the  plot  along  with  the  information 
listed  previously  with  FOURFIT  plots.  Figure  77  compares  the  data  record 
high  pass  filtered  at  the  fidelity  frequency  to  the  high  pass  equivalent 
Speicher-Brode  filtered  at  the  same  cutoff  frequency.  These  waveforms  are 
seen  to  compare  quite  well. 

The  results  of  high  pass  equivalency  have  not  been  adequately 
tested.  In  addition,  more  investigation  into  yield  dependence  is 


41 


warranted.  For  example,  Reference  10  discusses  a  method  for  removing  the 
phase  shift  resulting  from  a  filter.  Application  of  this  technique  may 
prove  useful.  Therefore,  though  the  work  is  promising,  the  FOURFIT  code 
as  presented  in  Appendix  A  does  not  include  the  capabilities  discussed  in 
this  section. 


SECTION  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

Ambiguity  and  uncertainty  in  performing  evaluations  of  airblast 
simulation  records  have  suggested  the  need  for  a  methodology  to  achieve 
consistent  and  meaningful  analysis  of  such  data.  Previous  studies  (Refs. 
1,  2,  3)  have  shown  that  the  FOURFIT  technique  meets  most  requirements. 

Until  the  present,  FOURFIT  has  been  used  to  graphically  determine  a 
best  fit  ideal  nuclear  waveform  for  a  simulation  record  based  upon 
comparison  of  the  data  Fourier  amplitude  spectrum  to  a  set  of  normalized 
Fourier  amplitude  spectra  derived  from  the  formulations  for  the  ideal 
nuclear  airblast  pressure  history  (i.e.,  "New  Brode"  or  Speicher-Brode) . 

In  addition  to  providing  consistent  estimates  for  equivalent  nuclear  yield 
and  peak  overpressure  for  records  from  a  single  event,  the  frequency 
analysis  provides  considerable  insight  into  the  frequency  content  of  the 
data  relative  to  the  ideal  nuclear.  Furthermore,  the  FOURFIT  methodology 
allows  for  determination  of  a  "fidelity"  frequency.  This  frequency 
indicates  a  cutoff  whereby  systems  with  frequency  sensitivity  at  or  below 
that  level  experience  a  good  simulated  loading  defined  by  the  equivalent 
nuclear  fit  determined  for  the  record.  Above  that  frequency,  the 
simulation  breaks  down. 

The  computer  code  FOURFIT,  and  its  companion  plotting  routine 
FOURPLT,  provide  an  automated  method  for  fitting  which  allows  for  rapid 
analysis  of  records  while  eliminating  analyst  bias.  In  addition,  the  code 
allows  for  studies  of  individual  records  and  of  Speicher-Brode  waveforms 
by  allowing  the  analyst  to  specify  a  fast  Fourier  transform  and  integrated 


impulse  or  low  pass  filtering  of  either  type  of  waveform.  These  codes 
were  written  for  use  on  the  Defense  Nuclear  Agency  CDC  Cyber  176  computer. 

The  results  of  the  applicaton  of  these  codes  for  analysis  of  two 
simulation  events,  0.35  KBAR  DISC  HEST  and  0.35  KBAR  HEST  indicates  that 
the  code  is  capable  of  determining  nuclear  fits  for  the  records  of  the 
former  of  these  events  which  compare  favorably  to  those  determined 
previously  using  the  graphical  methodology.  It  must  be  noted,  however, 
that  the  0.35  KBAR  DISC  HEST  data  base  consisted  of  high  fidelity, 

Spei cher-Brode-1 ike  pulses.  Data  from  the  second  event,  0.35  KBAR  HEST, 
were  not  of  such  high  fidelity.  These  records  were  analyzed  graphically 
in  the  study  of  Reference  3  and  at  the  time  were  found  to  be  difficult  to 
fit  due  to  their  poor  fidelity  marked  by  obvious  diversions  from  the  ideal 
nuclear  wave  shape.  The  FOURFIT  code  managed  to  fit  several  of  the 
records  rather  well.  However,  in  some  cases  the  poor  data  waveforms 
caused  the  fit  and  the  data  to  show  poor  agreement  at  late  time.  Closer 
examination  of  such  records  may  enable  better  fits  to  be  defined. 

However,  it  is  not  possible  to  automate  a  consistent  method  for  fitting 
non-typical  waveforms. 

Finally,  a  study  into  the  usefulness  of  high  pass  arid  band  pass 
filtering  yielded  mixed  conclusions.  Although  no  important  results  were 
recovered  from  the  band  pass  filter  study,  some  limited  insight  was 
provided  through  high  pass  filtering.  The  extent  of  this  effort  was 
limited  in  the  study  of  high  pass  filter  effects  and  so  was  not  totally 
conclusive.  However,  this  study  indicated  that  a  high  frequency 
equivalent  nuclear  waveform  may  be  estimated  through  application  of  high 
pass  filters. 
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Several  reconmendations  may  be  made  in  view  of  the  preceeding 
comments.  For  example,  when  the  available  data  from  an  event  proves  to  be 
of  low  fidelity  relative  to  the  ideal  nuclear  waveform,  a  means  for 
determining  guidelines  for  pursuing  the  fitting  of  such  data  needs  to  be 
developed.  Variable  weighting  schemes  for  the  squared  difference  values 
may  allow  the  analyst  to  better  address  the  different  power  regimes  in 
such  signals.  This  analysis  would  be  performed  on  a  case  by  case  basis. 

Next,  it  is  recommended  that  the  FOURFIT  code  be  applied  to  define 
record  fidelity  in  addition  to  the  low  pass  filter  definition  of  fidelity 
frequency.  It  is  suggested  that  the  methodology  may  be  expanded  to  enable 
quantification  of  fidelity  and  that,  with  increasing  interest  in  the 
development  of  a  high  fidelity  HEST,  a  set  of  fidelity  guidelines  may  be 
established  based  upon  a  scheme  of  sum  of  differences  between  the  data  and 
its  best  fit  Fourier  amplitude  spectrum.  Various  guidelines,  as  a 
function  of  frequency  range,  may  help  to  determine  the  relative  fidelity 
of  various  so-called  Hi-Fi  HEST  candidates. 

The  fits  to  the  normalized  Speicher-Brode  Fourier  amplitude  spectra 
have  only  been  checked  between  peak  pressure  values  of  1  MPa  and  200  MPa. 
There  is  increased  interest  in  higher  overpressure  regimes,  on  the  order 
of  600  MPa  and  above.  It  is,  therefore,  recommended  that  the  ability  to 
fit  simulated  overpressure  pulses  in  that  range  be  demonstrated  and/or 
developed.  This  would  require  a  study  of  the  Speicher-Brode  Fourier 
amplitude  equations  to  determine  additional  parametric  validity  up  to, 
say,  1000  MPa. 
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It  is  recommended  that  the  effects  of  high  pass  filtering  on  ideal 
and  simulated  nuclear  overpressure  histories  be  studied  further.  The  high 
pass  equivalency  technique  discussed  in  this  report  may  perhaps  be 
extended  to  evaluate  the  influence  of  high  pass  filters  on  estimates  of 
equivalent  yield  for  the  high  pass  fit.  This  might  be  accomplished 
through  use  of  the  "zero  phase  shift"  filter  as  discussed  in  Reference  10. 

Finally,  it  is  recoimended  that  the  computer  program  FOURFIT  be  used 
to  determine  equivalent  nuclear  yield  and  peak  overpressure  for  all  future 
simulation  events. 
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PE26.TAPE4e> 

PROGRAM  FOURFIT  ESTIMATES  THE  PEAK  OVERPRESSURE 
AND  NUCLEAR  YIELD  FOR  AIRELAST  SIMULATION  RECORDS 
BY  COMPARING  FITS  OF  THE  DATA  FOURIER  AMPLITUDE 
SPECTRUM  TO  THE  FOURIER  AMPLITUDE  SPECTRA  OF  TRIAL 
SPEICHER-6R0DES  RESULTS  ARE  WRITTEN  TO  A  FILE 
(  T  APE48  )  TO  EE  READ  AND  PLOTTED  E >  PROGRAM  FOURPLT 
WRITTEN  EY  C.  W  STEEDMAN.  APPLIED  RESEARCH  ASSOC. 
INC  .  A  L6UQUE  ROUE  .  NM ,  FEB  1984. 


COMMON  /FIlT 
COMMON  /PlCTV 


COMMON  /FFT  /  FRO( 3001 ). AMP( 3001 ) .XFFT ( 3001 ) 

COMMON  /'TERAT/  W< 5  )  . PI  5  )  .DELT AWl 5 ) . DELTAP! 5 ) . YLD! 5 ) 

COMMON  /THIST  /  TTIMI6000) . PRESSI6000) ,TIMP(2999 ), PIMP! 2999) . 
•  PF I L  T ( 6000 ) 

COMMON  /IMP  /  I  IMP .DTO , DTB . TPEB . DTBN 
COMMON  /POINTS/  NEPTS . NEPTS , N1 .NEF .NBF 
COMMON  /ESTIM  /  PSG1 . WI . PP . W13.PS0F , WF , F SO 
COMMON  /PEAK  /  DP . TA . PSO , ALPF 
COMMON  /S6C0NS /  RSKFT . YS . S . XM 

COMMON  /FIlT  /  I F I LT  . FLO( 7  I  . PFDM> ( 7  I . PFEM>  ( 7 ) 

COMMON  /PlCTV  '  1 T L ( t  1  . 1 S T L ( 6  )  . I DE 
COMMON  /UNITS  /  ’UNITS, JUNITS 
COMMON  /COUNT  /  1C0UNT . 10PT . LF ILT 
COMPLEX  XFFT 

t ape 2  contains  input  parameters 

NEPTS  *  NO.  OF  POINTS  TO  EE  READ  FROM  TAPE 
I  UN 1T$  -  1  FOR  TAPE  INPUT  PRESSURE  IN  PSJ 
*-1  FOR  TAPE  INPUT  PRESSURE  IN  MPA 
JUNITS  *  1  FOR  TAPE  INPUT  TIME  IN  MILLISECONDS 
■-1  FOR  TAPE  INPUT  IN  SECONDS 
PSOI  *  INITIAL  peak  OVERPRESSURE  ESTIMATE  IN  MPA 
WI  -  INITIAL  NUCLEAR  YIELD  ESTIMATE  IN  KT 
I  OPT  «  1:  FITTING  ROUTINE  TO  BE  DONE 

-  2:  JUST  FOURIER  TRANSFORM  THE  DATA 
«  3:  JUST  FOURIER  TRANSFORM  THE  SPE 1  CHE R -BRODE 
DEFINED  BY  PSOI ,  WI 
I F I L  T  *  1  FOR  FILTER  TO  BE  EXECUTED 
*  -  1  FOR  NO  F  I  LTER 

FLO  ■  LOW  PASS  CUTOFF  FREQUENCY  IN  HZ  (UP  TO  7  ALLOWED) 
(NOTE -FOR  LESS  THAN  FILTERS.  FLO 
MUST  BE  SET  TO  0.  TO  ESCAPE  THE  LOOP.) 


I  UN 1T$ 


REWIND  2 

READ( 2,111)  NEPTS, IUNITS. JUNITS 
RE  ADI  2 . 112)  PSOI ,WI 
RE  AD  (  2  ,  113)  IOPT.IFUT 
READ! 2.115)  ( FLO! I  )  , 1*1 , 7  ) 

111  F  ORMAT 1315) 

112  FORMAT! 2F5 . 2  ) 

1 13  FORMAT! 215  I 

1 15  FORMAT!  7F 10  0) 

WR I T  E i 6 ,  1  )  PSOI .WI 
1  FORMAT!  2X . -PSOI  «  *.F5.2.5X,*WI 
WR ITE I  48 .  113)  I  OPT , IF  I LT 
I  COUNT  •  0 
NBPTS  *  2048 
IF  I  I  OPT  £0.3)  GO  TO  7 
CALL  EBREAO 

IF! I0PT.E0.2)  00  TO  •«« 

CALL  FIT 
7  I  COUNT  •  1 
CALL  RANGE 
CALL  SPBROOE 
666  ENO 


* .F5.2  ) 


o  o  o  onn  non  n  nonnon  onooooo 


SUBROUTINE  EBREAD 


this  subroutine  reads  PRESSURE  VALUES  from  an 
EBCDIC  TAPE  BASED  UPON  THE  FORMAT  PREVIOUSLY 
USED  BY  WES. 


COMMON  /FFT  /  FRQf 300 1 ) . AMP ( 300 1 ) . XF FT ( 300 1 ) 

COMMON  /PGINTS/  NEPTS.NBPTS.NI.NEF.NEF 

COMMON  /THIST  /  TT I M ( 6000 )  . PRE SS ( 6000 ) . T I MP ( 2999  )  . P I MP ( 2999  ) 

•  PFILTI6000) 

COMMON  /FILT  /  I  F  I  L  T  .  F  LO  (  7  )  .  PFDMX  <  7  )  .  P  F  BMX  (  7  ) 

COMMON  /IMP  /  I  IMP . DTD . OTB . TPEB . OTBN 
COMMON  /UNITS  /  I  UN  ITS, JUN ITS 
COMMON  /PLOTV  /  ITl(8  I  .  ISTL(8  )  . IDB 
COMMON  /COUNT  /  I  COUNT ,  I  OPT . L F I lT 
COMPLEX  XFFT 

DIMENSION  IWKE ( 2000 ) . WKE ( 2000 ) 

DIMENSION  DUM ( 3  )  . DAI  5 ) 

DELP  IS  THE  DATA  BASELINE  SHIFT.  BE 
SURE  THAT  IT  is  IN  THE  PROPER  UNITS. 

DELP  *  0.0 
REWIND26 

READ  TAPE  HEADER  INFORMATION 

RE  AD ( 26 , 30 )  I T  L ( 3  )  . ITL(4  )  . 

•  DUM (  1  ) . DUM ( 2  )  . 

•  I T  L  C  1  ) .  I TL  < 2  )  . 

•  DTD.NP 

30  FORMAT! 3! 2A 10  )  . E 15  e . IS t 

ITUS!  *  1 0H  PRESSJRE 
I  T  L.  (  6  )  -  10HHISTORY 

I T  L ( 7  )  «  1 0H 

I T  L ( 8  )  *  1  OH 

WRITEI4S.35)  ( I TL < L ) . L *  1 . 8  ) 

35  FORMAT ( BA  10  ) 

DO  20  1*1 .NEPTS 
TT I M(  I  )  =  0. 

PRESS! I  )  *  O. 

20  CONTINUE 

I F ( E OF ( 26  )  )  900.901 

SET  UP  DATA  UNITS  CONVERSIONS; 

MSEC  TO  SEC  AND  PSI  TO  MPA. 

SOI  IFIJUMTS.GE.  1)  DTD  *  DTO-.OOI 
PFACT  *  .006894757 
IF ( IUNITS . LT .0)  PFACT  •  1. 

IP  *  1 

time  *  o. 

NLINE  -  NEPTS/5 
RUNE  *  FLOATfNEPTS  >/5. 

IF ( RLINE . GT . NLINE )  NLINE  •  NL I NE ♦ 1 

READ  PRESSURE  VALUES 

DO  40  J*1 .NLINE 

READ! 26.50)  ( DA< JJ ) , JJ* 1 . 5  ) 

50  F0RMAT(5E 16.8) 

I F ( EOF ( 26  i  )  900.902 
•02  DO  60  K* 1 . 5 

TTIM! IP |  *  TIME 
P  *  DA ( K  ) 

PRE S S I  IP  !  *  ( P-PFACT  ) -DELF 
IP  =  I P-  1 
TIME  *  TIME-DTD 
60  CONTINUE 

40  CONTINUE 


SPLINE  THE  END  OF  THE  DATA  TO  ZERO  IN 
CASE  Or  A  TRUNCATED  RECORD 
TLAST  .  TTIMiNEFTSj 

CALL  SPL INE ( TLAST . NEPTS . TTIM . PRESS  I 
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C  IF  JOPT  •  1,  FIND  THE  TIME  TO  DATA 

C  PEAK  TO  AID  IN  PHASING  THE  OVERLAYS. 

c  aid  in  phasing  overlays 

PMAX  *  0. 

DO  78  I K  *  1  , NE PT S 

PMAX  *  AMAX 1 ( PMAX . PRESS! IK )) 

I F ( PM A  X  EQ. PRESS! IK)  )  T  PE  B  *  TTIM(IK) 

78  CONTINUE 
C 

c 

C  REMOVE  BASELINE  CORRECTION  FOR  POINTS 

C  BE FORE  the  ARRIVAL  OF  THE  SHOCK 

DO  77  M= 1 .NEPTS 

IF(TTIMiM)  GT  TPEEI  GO  TD  990 
PRESSIMI  =  PRESS! M l+DELF 
77  CONTINUE 
C 

GO  TO  990 
900  WRj  TE ( 6 . 70  > 

70  FORMAT ! 10* , * END -OF -FILE  REACHED  EARLY*..,/) 

990  CONTINUE 

I F (  I  OPT . NE  . 2  >  GO  TO  45 

CALL  F MAXIPRESS. NEPTS. YPMN, YPMX ) 

CALL  FMAX! TT IM . NEPTS . XPMN. XPMX ) 

WRITE ( 4e . 100)  NEPTS . XPMN . XPMX . YPMN . YPMX 
IF ( IFILT . LT  0)  GO  TO  700 
C 

C  CALL  FOR  FILTERS  TO  BE  EXECUTED 

CALL  FLOOPfTTIM.PRESS.DTD.NEPTS.PFILT  ) 

RETURN 

C 

700  WRITEI48. 105)  (TTIM(K).K*1, NEPTS) 

WRITE ( 4e .  1 QS )  ( PRESS! KL  )  ,KL= 1  .NEPTS ) 

100  FORMAT!  15 . 4E 15 . 8  ) 

105  F0RMAT ( 1 0E  '  5  8  ) 

4  5  I  IMP  «  1 
C 

C  IMPULSE 

C 

CALL  IMPULSE!  IIMP.OTD. NEPTS. NI  ) 

I F  c I  OR  T .NE .2)  GO  TO  110 
1  T  L  I  5  t  *  10H  IMPULSE  H 
I  T  L(  6  )  *  1 0H  I  ST  OR  Y 
WR ITE(48,  115)  ITL(5  )  .  I TL ( 6 ) 

CALL  FMAXITIMP.NI .XIMN.XIMX ) 

CALL  FMAX.IPIMP.NI  .Y1MN.YIMX) 

WRITE (48 .  100)  NI  . X  I MN . X  I  MX . Y I MN . Y I  MX 
WRITE ( 48 ,  105  I  (TIMPI  1H).IH*1  .NI  ) 

WRITE ( 4e . 105 )  < PIMP! JH ) . JH* 1 . NI ) 

115  F0RMAT(2A10> 

C 

C  find  the  FOURIER  TRANSFORM  AND  CALCULATE  AMPLITUDE. 

C 

1 10  TTOT  *  DTD* NEPTS 
C  FREQUENCY  INCREMENT 

DFE  *  1  .  /TTOT 

FOE  *  0. 

C  fqurieR  TRANSFORM 

CALL  FFTRC! PRESS. NEPTS. XF FT, I WKE.WKE ) 

XRE  •  REAL! XFFT( 1 ) )/(2*NEPTS ) 

X] E  •  AIMAGIXFFT! 1 ) )/(2*NEPTS) 

FOE  •  FQE*DFE 
C  AMPLITUDE  SPECTRUM 

FRO! 1 )  ■  FOE 

AMP ( 1 )  *  SORT! 2 . • l XRE • XRE*X1E*XI E ) ) -TTOT 

NEF  *  NE PT S  ' 2*  1 

DO  eo  UK* 2 . NEF 

c  OE  =  FQE*DFE 
PRQI JK I  *  FOE 

>RE  *  REAL! XFFTI JK ) l/NEPTS 
XIE  *  AIMAG!  xFFT! JK  )  i/NEPTS 

AMP(JK)  •  SORT(XRE'XRE*XIE*XIE )*TTOT 
80  CONTINUE 
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IF ( IOPT.NE . 2)  RETURN 
1TLI5)  •  10H  FOURIER  A 
I T l ( 6 )  *  10HMPLITUDE  S 
ITL(7>  =  10HPECTRUM 
CALL  FMAX ( FRO . NEF . XFMN .XFMX ) 

CALL  FMAX( AMP . NEF . YFMN, YFMX ) 

WR I T  E ( 4  B ,  117)  ITL(5),ITL(6).ITL(7) 

117  F  ORMAT ( 3A 10 ) 

WR I T E  (  48 .  lOO )  NEF . XFMN . XFMX . YFMN, YFMX 
WRITE(48.  105)  ( FRO( L I  )  . L I ■ 1 . NEF ) 
WRITE(48. 105)  ( AMP ( JI ) . JI ■ 1 . NE F ) 
RETURN 
END 

SUBROUTINE  FIT 


THIS  SUBROUTINE  ITERATES  ON  YIELD  WITHIN  ITERATIONS  ON 
PEAK  PRESSURE.  ITS  AIM  IS  TO  REDUCE  THE  SUM  OF  THE  SOUARES 
OF  THE  DIFFERENCE  BETWEEN  THE  DATA  AMPLITUDE  AT  F ( I )  ANO 
THE  ESTIMATED  SPE1CHER -8 RODE  AMPLITUDE  AT  F(I)  DIVIDED 
BY  Ft  I)  BASED  UPON  A  TOLERANCE  ON  PEAK  PRESSURE  AND  YIELD. 
END  RESULT  IS  A  FINAL  ESTIMATE  OF  PEAK  OVERPRESSURE 
( PSOF  )  AND  YIELD  (WF).  ALSO,  AN  ESTIMATE  OF  THE  GOODNESS 
OF  FIT  (DELL)  IS  DETERMINED.  PRESSURE  IS  IN  MPA,  YIELD  IS 
IN  KT  . 


COMMON  /POINTS/ 
COMMON  /ESTIM  / 
COMMON  /ITERAT/ 
COMMON  / F  F  T  / 
COMMON  /PEAK  / 
DATA  TOL/.OI/ 


NEPTS.NBPTS.NI .NEF ,NBF 

PSOI , WI , PP , W13 , PSOF . WF . FSO 

W<  5 ) , P ( 5 ) , DELTAWI 5 ) , DE L T AP ( 5  )  , Y LD ( 5  ) 

FRQ( 3001 ) .AMP (3001 ) . XFFT( 3001 ) 

DP.TA.PSO.ALPF 


P( 1 )  *  . 1-PSOI 
P ( 2 )  =  ,4-PSOI 
P(3)  =  1.0-PS0I 
P ( 4  )  *  4. -PSOI 
P ( 5  )  =  10.-PS0I 
JPRESS  -  0 


LOOP  ON  PRESSURE  TOLERANCE 

DO  100  UU  = 1 .50 

UPRESS  -  JPRES'.+  I 
UMIN  *  2 
JMAX  =  4 

I F ( UPRESS . NE  1  )  GO  TO  105 
UMIN  «  1 

UMAX  -  5 


LOOP  ON  PRESSURE 

105  DO  200  I  I -UMIN, UMAX 
PP  •  P (  I  I  ) 

UYLO  *  0 
W( 1 )  *  0.1 *WI 
W ( 2 )  *  0. 4-WI 
W<  3  I  «  1 . 0*WI 
WI  4 )  =  4 . 0*WI 

W(  5 )  ■  10 . *WI 

LOOP  ON  YIELD  TOLERNACE 

DO  250  KK* 1 . 50 

U  y  LD  *  UYLD+’ 

IMIN  ■  2 
IMAX  •  4 

IF ( UYLD.NE . 1 )  GO  TO  255 
IMIN  ■  1 
IMAX  •  5 

LOOP  ON  YIELD 

555  DO  300  LL* IMIN , IMAX 

W 1 3  •  W(  LL  )  • ■  . 333C3 
1F(LL  NED  GO  TO  256 
CALL  RANGE 


52 


WJL 


*  */  m f  * _•  »  * _  ■| 


oou  ouu  uou  o u u o ooooooo u o 


determination  of  residuals 

256  DELTAW) LL )  *  0. 

DO  350  LK« 1 .NEF 

F  SC  L  *  FRQ(LK)‘W13 
I F ( FRQ( LK) .GT . 7000. )  GO  TO  300 
I F I F SCL  LT . F SO )  GO  TO  350 
CALL  AMPALGt  FSCL , BAMP ) 

AMPN  *  ALOG 10( AMP ( LK ) ) 

BAMPN  •  ALOGIO(BAMP) 

DF2  •  FRO (LK)*FRQ(LK) 

DELTAA  *  ( AMPN - BAMPN  > /FRO( LK) 

IF ( FRQI LK ) . LT . 1000. )  DELTAA  *  2. ‘DELTAA 
IF ( FRQI LK ) . GT . 5000.  .AND.  FRQ( LK ) . LT . 7000 . ) 

•  DELTAA  *  2. ‘DELTAA 

DELTAA  «  DELTAA‘D£LTAA 
DELTAW(LL)  ‘  DELTAW(LL)+DELTAA 
350  CONTINUE 
300  CONTINUE 

RESET  YIELDS 

EPSW  *  AES( WI 5  )-W(  1 ) )-2 . /( WI 5 l»WI  1  )  ) 

I F ( EPSW . LT  TOL )  GO  TO  360 
CALL  RESETW 
250  CONTINUE 

WR I T  E I  6  .  1250> 

1250  FORMAT ( 2X . *FAI LED  TO  CONVERGE  ON  YIELD") 

STOP  14 
360  CONTINUE 

DWMIN  *  AMIN1 (DELTAWI  1  )  .OELTAWI 2 ! .DELTAWI 3 ) .OELTAWI  4 ) . DELTAWI  5 ) ) 
DO  365  MM" 1 ,5 

IFIDELTAWIMM) . £Q. DWMIN)  KW  *  MM 
365  CONTINUE 

Y  LD (  II)  *  WIKW  ) 

DELTAPI II)*  DELTAWIKW) 

200  CONTINUE 

RESET  PRESSURES 

EPSP  *  ABS< P( 5 l-P<  1  )  )‘2 - /( P<5  )*P(  1  )  ) 

1FIEPSP  LT  TOL)  GO  TO  400 
CALL  RESETP 
100  CONTINUE 

WRITE (6. 1100) 

IIOO  FORMAT (2X.  ‘FAILED  TO  CONVERGE  ON  PEAK  PRESSURE •  ) 

STOP  10 

400  DPMIN  •  AMIN1 ( DELTAPI 1  )  .DELTAPI 2  )  .DELTAPI 3 ) .DELTAPI 4 ) .DELTAPI 5 ) ) 
DO  405  NN* 1 . 5 

IFIDELTAPI NN) . EO. DPMIN )  KP  *  NN 
405  CONTINUE 

W 1  3  *  YLD(KP)“  .33333 
pp  .  P(KP ) 

DELL  *  DELTAPIKP )/NEF 

RETURN 

END 

SUBROUTINE  AMPALGI FSCL .BAMP ) 


THIS  SUBROUTINE  ESTIMATES  THE  FOURIER  AMPLITUDE  OF  THE  TRIAL 
PEAK  PRESSURE  AND  YIELD  BASED  UPON  A  FIT  TO  THE  SUITE  OF 
OF  NORMALIZED  SPEICHER-BRODE  FOURIER  AMPLITUDE  SPECTRA. 

THE  ALGORITHM  USES  SCALED  FREQUENCY  OF  INTEREST  (FSCL). 
SCALED  FUNDAMENTAL  FREQUENCY  OF  THE  S'6  OF  CONCERN  I F SO ) 

AND  THE  PEAK  OVERPRESSURE  (PP)  TO  CALCULATE  THE  SCALED 
AMPLITUDE.  THE  ALGORITHM  USES  PRESSURE  IN  MPA  AND  YIELD 
IN  KT .  THE  EQUATIONS  ARE  FOR  A  SURFACE  BURST  ONLY.  THEY  ARE 
VALID  FOR  ANY  YIELD  AND  FOR  PEAK  OVERPRESSURE  UP  TO  10OMPA 


COMMON  /ESTIM/  PS01 ,¥I , PR , W13 . PSOF , WF , F50 

At  •  .  1 788‘PP*  * ( - . 72 ) • <  FSCL  “  ( -  1 . »PP» • (-.103)1) 

A2  •  ,01474‘PP“(  -  .  15  )•(  FSCL/FSO  )“<  -  1  .  75  ) 

A3  •  .  001  1  *PP“  I  PP‘ •(  -  .  234  i  )•(  FSCL/FSO  )  -  •  I  -2  .  15  ) 

A4  *  .  00132"FSCL“(  -  .  547  ) 

AS  •  . 01034-PP* • I - .  1  13  )• I  1  .  'FSCL I • I FSCL/FSO  1“  I  -  1 . 5  I 
A6  •  .OOOOl  1  *PP“  .  77‘ (  FSCL/FSO  )  “  I -7 . 5  ) 

A 7  *  ,0000€66‘PP“  .  3‘(FSCL/FS0)“(  -  1 . 5) 

ASCL  •  A  1  -42*43  +  44  +  45-46*47 

BAMP  •  ASCL‘PP‘W 1 3 

RETURN 

END 


nnnn  no  oooooo  no  nn  on  no  noon  on  nnonoo 


subroutine  reset* 


THIS  SUBROUTINE  RESETS  THE  FIVE  YIELD  VALUES  BASED 
UPON  THIS  ITERATION'S  MINIMUM  RESIDUAL 


COMMON  /ITERAT/  W< 5  )  . P < 5  )  . OE L T AW) 5  )  . DE L T AP < 5  )  . Y LD < 5  ) 

FIND  THE  MINIMUM  DELTA 
IF(DELTAW) 5) . LT . DELTA*) 4 )  )  GO  TO  10 
IF(DELTAW(4 ) . LT . DELTA*) 3 ) )  GO  TO  20 
IF(DELTAW) 3 ) . LT .DELTA*) 2 ) )  GO  TO  30 
IF(D£LTAW( 2 ) . LT .DELTA*)  1 )  )  GO  TO  40 

REDEFINE  YIELDS  BASED  UPON  THE  MINIMUM 

IF  DELTA*) 1 )  IS  MIN. 

DYLD  =  (W(2t-*<  1 )  )  *0 . 25 
*(  5  )  *  W(  2  ) 

DELTA*) 5 )  =  DELTA*) 2  ) 

GO  TO  50 

IF  DELTAW(5)  IS  THE  MINIMUM, 

10  DYLD  =  ) W( 5  I -*( 4 ) ) -0. 25 
W(  1  )  =  W(  4  ) 

DELTA*)  1  )  =  DELT  A*) 4 ) 

GO  TO  50 

IF  DE  LT  AW ( 4 )  IS  THE  MINIMUM, 

20  DYLD  =  (W(5 )-*( 3) )»0. 25 
*(  1  )  *  W( 3  ) 

DELTA*)  1  )  -  DE  LT  A* ( 3 ) 

GO  TO  50 

IF  DELTAW<3)  IS  THE  MINIMUM. 

30  OYLO  *  (W(4)-*)2) >*0.25 
*11)  *  W( 2  ) 

*15)  =  ta( 4  ) 

OE  ltawi  1  )  =  DELTA*) 2  I 
DE  l.  T  A  w  i  5  )  -  DELTAW)  4  ) 

GC  TO  50 

IF  DElTAW(21  IS  THE  MINIMUM, 

40  D>  LD  =  ( W< 3  )-*(  1 )  )*0. 25 
W(  5 )  =  *( 3  ) 

DELTA*(5)  *  DELTA*) 3  ) 

50  *( 2 )  =  *(  1  )*DYLD 
*( 3  )  *  w( 2  )*DYLD 
*( 4  )  •  *13  )*DYLD 
RETURN 
END 

SUBROUTINE  RESETP 


THIS  SUBROUTINE  RESETS  THE  FIVE  PRESSURE  VALUES 
EASED  UPON  THIS  ITERATION'S  MINIMUM  RESIDUAL. 


COMMON  /ITERAT/  W( 5  )  . P < 5  )  . DE LT AW ( 5  )  . DE LT AP ( 5  )  , Y LD ( 5  ) 

find  the  minimum  deltap 

IF) CELTAP) 5 ) . LT . DELTAP) 4  )  )  GO  TO  10 
IF ( DELTAP) 4 )  . LT .DELTAP) 3  )  )  GO  TO  20 
IF) CELTAP) 3  )  . LT ,DELTAP(2  )  )  GO  TO  30 
IF)CELTAP(2). LT. DELTAP)  1  ))  GO  TO  40 

REDEFINE  PRESSURES  BASED  UPON  THE  MINIMUM 


,V, 

> 


IF  DELTAP) II  IS  THE  MINIMUM. 
DPRESS  ■  ( P ( 2  )  *  P (  1  I  )  *0 ■ 25 
P(  5  )  •  P<  2  ) 

*(  5  )  ■  *(  2  I 

DE  L  T  AP ( 5  I  *  DELTAP) 2  ) 

GO  TO  50 


54 


c 

C  IF  DELTAP(E)  IS  THE  MINIMUM. 

10  DPRESS  *  (P<  5)-P<4  )  )•  . 25 

P (  1  )  *  P(4  ) 

W(  1  )  *  W(  4  ) 

DElTAP(  1  )  *  D£LTAP(4  ) 

GO  TO  50 
C 

C  IF  DE IT AP ( 4  )  IS  THE  MINIMUM. 

20  DPRESS  *  ( P ( 5 ) - P ( 3 ) ) «0 . 25 
P< 1 )  *  P<3) 

W(  1  )  «  W(  3  ) 

DELTAPf  1  )  *  DE LTAP ( 3  ) 

GO  TO  50 
C 

C  IF  DELTAPO)  IS  THE  MINIMUM. 

30  DPRESS  *  (P(4)-P(2))»0.25 
P  (  1  )  *  P  (  2  ) 

W  (  1  )  *  W  (  2  ) 

DELTAPl 1 )  *  DE  LT  AP ( 2 ) 

PI  5  I  *  P  (  4  ) 

W(  5  )  ■=  W(  4  ) 

DE  L  T  AP ( 5  I  =  DELTAP14) 

GO  TO  50 
C 

C  IF  DE  L  T  AP ( 2  )  IS  THE  MINIMUM. 

40  DPRESS  =  (P(3)-P( 1 ) )*0.25 
P(  5 )  =  P<  3  ) 

W<  5 )  =  W( 3  ) 

DE  L  T  AP ( 5  )  *  DE  L  T  AP ( 3  ) 

50  P( 2 )  =  P(  1  l  +  DPRESS 
P ( 3 )  *  P ( 2  t+DPRE  SS 
P(4)  =  P ( 3 (+DPRESS 
RETURN 
END 

SUBROUTINE  RANGE 


C  THIS  SUBROUTINE  IS  AN  ITERATION  TD  FIND  THE  RANGE 

C  OF  THE  ESTIMATED  PEAK  PRE5SURE  FOR  THE  ESTIMATED 

C  YIELD.  THIS  IS  NECESSARY  FOR  COMPUTATION  OF  THE 

C  SPEICHER-BRODE  PRESSURE  HISTORY.  TIME  OF  ARRIVAL 


on  o  non  non  o  ooooooo 
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SUBROUTINE  PP£AK( K . V . PE AKP  ) 


THIS  SUBROUTINE  CALCULATES  THE  PEAK  OVERPRESSURE  (MPA), 
TIME  OP  ARRIVAL  ( T A , MS/KT • •)/?).  AND  POSITIVE  PHASE 
DURATION  ( DP , MS/KT  * • 1 /3  )  AFTER  SPE I  CHER -BRODE ,  JUNE.  1982. 


COMMON  /PEAK  /  DP . TA . PSO . ALPF 
COMMON  /SbCONS/  RSKFT.VS.S.XM 

XLEAST  •  1.E-9 
VLEAST  *  1.E-9 
2MA  X  *  100. 

I F I  A . LT . XLEAST  )  X  ■  XLEAST 
IF( Y . LT . VLEAST  )  Y  «  YLEAST 


R  « 

SORT  (  x 

R2  * 

R*R 

R3  * 

R*R2 

R4  * 

R2-R2 

R6  * 

R  2  •  R4 

R8  • 

R4  •  R4 

2  * 

Y/X 

22  = 

2*2 

23  * 

2-22 

25  ‘ 

22-23 

217  *  2* •  17  . 

218  *  2*-  18. 

V7  -  Y*  -7  . 

IF (2  .  GT . 2MAX  )  2  -  2MAX 

XM  *  170. • Y/{ 1 . +337 . • Y* • . 25 >•* . 914- Y**2 . 5 

SCALED  TIME  OF  ARRIVAL 

U 1  *  (  .543-21  .B-R+386. *R2+2383. *R3)*R8 

U2  *  2  99E- 14- 1 . 91E - 10-R2* 1 .032E-6-R4-4 . 43E-6-R6 

U3  *  ( 1 .028*2  087- R+2.69-R2 ) • R8 

UTA  *  U  1  / ( U2  +  U3 ) 

TA  *  UTA 

IF ( X . LT . XM)  GO  TO  101 

W1  «  < 1 .086-34  605-R+486  3-R2+23B3. *R3 >*R8 
W2  *  3.0137E- 13- 1 . 2 128E-9- R2+4 . 128E-6-R4- 1 . 1 16E -5-R6 
W3  *  (  1 .632  +  2  629-R+2  69-R2  1-Re 
WT A  *  W  1/(  W  2  +  W  2  ) 

TA  *  UTA-XM/X+WT  A* (  1 . -XM/X  ) 

SCALED  POSITIVE  PHASE  DURATION 

S  *  1  .  -  1  .  IE  10* Y7/(  1  .  *1 .  IE  10* Y7 ) - ( 2 . 44  IE -8* Y* V/ 
(1.+9.E10*Y7))*(1./(4.41E-11+X**10. )) 

DP  *  <(  1640700. *24629  *TA  +  4  16 .  15-TA-TA )/ 

( 10880. *6 19 . 76*TA+TA*TA) ) 

*<  . 4+  . 001204* (T A* • 1 .51/(1*. 00 1 559* T A** 1.5)+ 

•  (  .04  26+ .5486*(TA**.25)/(1.  +  . 00357*TA*  *  1 . 5 ) )*S ) 

AA  •  1 . 22- (3. 908*22 )/( 1 . *810. 2*25) 

BB  *  2.321*<Zl8/(  1  *1.  113-Tig)  )*©.  195-(  . 038  3 1*217)/ 

(  1  + . 024  15-217)* . 6692 /( i . +4 164 . *2*  *8 .  ) 

CC  *  4 .  l53-<  1  .  149-218  )/(  1  ♦  1  .641-218)- 1 .  1  Ml  .  +  2.771*2.-3.5 i 
DD  «  -4 .  166* ( 25. 76- 2*' 1 . 75  I  M . ♦ 1 . 3e2- 2 18  1  +  6 . 257*2/  (  1  .  +  3 . 2  19*2  I 
EE  --  1  -  (  004  64  2  -2181/(1.*  003886-  Z  1 E  i 

FF  «  . 6096*(2. 879*2**9.25)/ ( 1  *2 . 359*2** 14 .5)- 17 . 15*22/ 

(  1. +71. 66-23) 

GG  •  1 . 83*5 . 361 *22/(  1 .♦. 3139*2* *6 .  ' 

HH  •  - (64 .67*25+ . 2905 )/(  1  *44  1 . 5*25 '-  1 . 389*2/(  1 . +49 .03-25)  + 

( 8 . 808- 2*  *  1 . 5  >/(  1 .  +  154  5*2*  *3  5 )+(  . 0014*R2/(  1 . -  .  158-R  + 

. 04  86-fi* •  1  5+ . 00128»R2 ) )*( 1 . /(  1  . +2 . • Y  )  ) 

PEAK  OVERPRESSURE 

PO  •  10.47/(R.»AA )+BB/(R**CC )+DD»EE/( 1 . +F F • R» « GG >+HH 

PEAKP  »  P0+ .00***4757 

RETURN 

END 
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SUBROUTINE  SPBRODE 


•  • 


C  THIS  SUBROUTINE  calculates  the  pressure  history  for 
C  THE  FINAL  PRESSURE-YIELD  PAIR  DETERMINED  BY  SUBROUTINE 

C  FIT.  IT  USES  THE  S°EICHER-BRODE  JUNE. 1982  ALGORITHM. 


COMMON  /THIST  / 

COMMON  /FFT  / 
COMMON  /ESTIM  / 
COMMON  /PEAK  / 
COMMON  /FILT  / 
COMMON  /SBCONS/ 
COMMON  /POINTS/ 
COMMON  /IMP  / 
COMMON  /COUNT  / 
COMMON  /PLOTV  / 
COMPLEX  XFFT 


TTIM(  6000 ) , PRESS ( 6000 ).TIMP (2999), PIMP (2999) 
PF I LT ( 6000 ) 

FRO( 3001 ) .AMP ( 3001 ) .XFFT (3001 ) 

PSOI ,W1 . PP.W13 . PSDF ,WF . FSO 
DP.TA.PSO, ALPF 

IFILT .FLO (7 ) ,PFDMX(7) .PFBMX(7 ) 

RSKFT.YS.S.XM 

NEPTS . NBPTS . NI .NEF ,NBF 

IIMP.DTD.DTB.TPEB.DTBN 

ICOUNT , IOPT , LFILT 

ITL(8).1STL(8). I OB 


DIMENSION  I WKE ( 1 1 ) 
DATA  JCOUNT / 0/ 


IF< IOPT .NE . 3 )  GO  TO  5 
ITl(1)  -  lOHCALCULATED 
I  T l ( 2  I  =  1 0H  SPEICHER- 
I T l ( 3  )  *  10HBRODE  PRES 
ITU  4  )  =  10HSURE  HI  STO 
I  T  L ( 5  )  *  lOHRY 
I TL ( 6  )  -  10H 

I T  L ( 7  )  *  10H 

I TL ( 8  )  *  10H 

WR I TE ( 48 , 26  )  (  I  TL  (  I  0 ) . 1 0» 1 . 8  ) 
26  F  ORMAT ( 8A 10 ) 


200 

210 


CALCULATE  SPE I  CHER -BRODE  TIMESTEP  BASED 
UPON  THE  POSITIVE  PHASE  DURATION. 

DTB  »  DP/NBPTS 
GO  TO  15 

ISTL( 1 )  «  lOHWITH  FOURF 
I  STL ( 2 )  •  lOHIT  SPE I  CHE 
I  STL ( 3 )  *  10HR-BROOE 
ISTL(A)  *  1 OH 
I  STL ( 5  )  *  10H 
I  ST l  (  6  )  *  1 0H 
I  ST  L ( 7 )  *  1 0H 
I  ST  L ( 8 )  -  1 0H 

WR I TE ( 48 . 26 )  ( ISTL( IG) . IG* 1 .8) 

CALL  FMAXI PRESS . NEPTS . YPMN. YPMX ) 

CALL  FMAX< TTIM, NEPTS . XPMN . XPMX ) 

WRI TE ( 48 . 200 )  NEPTS . XPMN , XPMX , YPMN . YPMX 
WR I T E ( 48 . 2 10 )  ( TTIM( IU ) , IU* 1 .NEPTS  ) 

WRITE (48. 210)  (PRESSt IP  )  , IP* 1 .NEPTS) 
FORMAT ( 15. 4E 15.8 ) 

F  ORMAT ( 10E15.8) 

ICOUNT  *  O 


FIND  THE  PEAKS  OF  THE  LOW  PASS 
FILTERED  DATA  PRESSURE  HISTORIES 
DO  7  1*1,7 

CALL  FILTER(DTD. NEPTS) 

CALL  FMAX(PFILT, NEPTS, PFDMN , PFDMX ( I ) ) 
7  CONTINUE 
ICOUNT  *  1 


CALCULATE  SPE I  CHER  -BROOE  TIME  STEP  BASED 
UPON  THE  DATA  TIME  STEP  FOR  FILTERING 
DT&  *  DTD- 1000. /W 13 

CALCULATE  THE  SPEICHER-BROOE  TIME  STEP  BASED 
UPON  THE  POSITIVE  PHASE  DURATION  FOR  OVERLAYS 
35  I F ( JCOUNT . EO . 1 )  DTB  •  DP/NBPTS 
15  DO  25  KJ* 1 , NBPTS 

TTIM(KJ)  •  0. 

PRESS(KJ)  *  0. 

25  CONTINUE 


'.r  ip* •  ».>  j  1  1  1  r  <  i.  ■  .-  ->  • 


C 

c 

c 


c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


X  «  R$KFT 

T  F  ■  TA*OP 

PO  •  PSOF • 145 . 038 

F  «  (  .01477- ( TA»  * . 75 )/( 1 .  +  .005836. T A )  +  7 . 402E - 5* ( T A  * • 2 . 5)/ 

•  <  1  *1  429E -8*TA**4  75 >- . 2 16  I *S+ . 7076-3 .077E-5* 

•  TA*TA*TA/( 1 ,+4.367E-5*TA*TA*TA) 

G  *  10 .+( 77  58-64 . 99* ( TA** .  125 )/( 1 .  +  . 04 348 • SORT C T A  )  )  )*S 
H  «  2.753*.05601*TA/M .  +  1 . 4 7 3E - 9 • T A • *5 .  )■*(  ,01769*TA/ 

•  ( 1 . *3 . 207E - 10*TA**4 . 25 ) - . 03209* (TA*»1 ,25)/( 1 .♦9.914E-8* 

•  T A* • 4 .  >-  1 . 6  )*S 

CALCULATE  PRESSURE  HISTORY 


DO  400  0*1 .NePTS 

T  *  TA+(J-2)*DTB 
SAVE  UNSCALED  TIMES 

TTIM(J)  ■  T  *W 1 3/ 1000 . 

PRESS! 0 )  ■  O. 

if<t.lt.ta)  go  to  aoo 

I  F  <  T  GT , TF  )  GO  TO  410 

B  •  ( F • ( TA/T  )  - - G*  <  1  -F  I* ( TA/T )**H)*(  1  . -( T-TA )/DP  ) 

POFT  *  PO-E 

JFO.LT.XM.  OR  .  Y  GT  0.38)  GD  TO  390 
XE  *  3  039- 1/(1. ♦6  7- y  | 

E  *  AES  I  l >  -  XU  1  / I >  E -XM  (  I 
IF ( E - GT . 50.  )  E  *  50. 

D  *  .  23-*5e3000  26667  .♦  1  .  E6*Y*  Y  )•*.  27*  E  + (  .  5-583000  V  •  Y/ 

•  ( 26667 . * i . E6» y » y )  )*E • *5 . 

A  *  (  r  -  1 .  1 • ( 1 . - ( E ■ ■ 20 .  I  / (  1 . ♦ ( E  *  *  20  ) ) ) 

07  »  474 . 2*  Y  *0  -XM  •  1 . 25 
I F ( DT . lT  .  1  .  E -9  )  DT  *  1  . E-9 
GA  =  (T-TA)/D7 
I F ( GA . GT , 400 . )  GA  *  400. 

v  *  1 . ♦( 3 . 28E 1 1 • ( Y»*6 . )/( 1 . *  1 . 5E 12* Y» *6  75 ) )* ( GA*GA*GA/ 

.  (6 . 13+GA-GA-GA ))*<1./(1.+9.23*E*E>) 

C  *  ((  1  04-240. 9* ( X* *4 )/(  1 . +231 . 7*X* *4  )  )• ( GA* *7  )/ 

.  ((  1.*.923*GA**85)*(  1  ♦A)  ) )*( 1 . - ( ( T-TA  )/DP  )**8  .  ) 

•  *2 . 3E 13* Y**9. /( 1 . *2 . 3E 13* Y--9  ) 

POFT  *  P0*( 1  +A  )• (B*V*C  ) 

390  PRESS(J)  *  POFT /  145 . 

400  CONTINUE 

410  UCOUNT  »  JCOUNT+1 

UNSCALE  THE  SPEICHER-BRODE  TIMESTEP 
DTBN  *  DTB-W13/ lOOO. 

IF ( JLOUNT . GT . 1  .OR.  I0PT.E0.3)  GO  TO  900 

FIND  THE  PEAKS  OF  THE  LOW  PASS  FILTERED 
SPEICHER-6R0DE  PRESSURE  HISTORIES 
LFILT  *  0 
DO  17  0*1,7 

CALL  FlLTERfDTEN. NEPTS) 

CAVL  FMAX ( PF I lT .NBPTS . PFEMN. PFBM/  (  J )> 

17  CONTINUE 

FIND  THE  LOW  PASS  FIDELITY  FREQUENCY 

DO  27  K* 1 . 7 

PFMAX  *  PFDMX(K(*0  90 
I F ( PFMAX .LE.PFBMX(K))  GO  TO  47 
27  CONTINUE 

WR2  TE ( 6 , 37  ) 

37  FORMAT ( 2X , •  +  ♦+  FAILED  TO  LOCATE  LOW  PASS  FIDELITY  ) 

ALPF  *  -999. 

WR I T E ( 4 8 . 57 )  ALPF 
GO  TO  35 

47  ALPF  *  FlO(K) 

WRITE  (  48 .5*7  )  ALPF 
57  FORMAT* F 10. 0 ) 

WR1  TE ( 6 , 67  )  ALPF 

67  FORMAT ( 2X, LOW  PASS  FIDELITY  (HZ)  •  *.F10.0.*  ♦*♦•) 

IF ( UCOUNT . EO. 1 )  60  TO  35 


59 
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DETERMINE  NUMBER  OF  SPEICHCR  BRODE  PAIRS  TO 
BE  PLOTTED  FOR  OVERLAY 
900  TE  •  NEPTS'OTD 

NPPTS  «  IFIX(T£/DTBNI 
1 F ( I  OPT . EG  -  3  )  NPPTS  *  NBPTS 
MR  I TE ( 48 , 4  50 )  NPPTS 
450  FORMAT! 15) 

I F ( IOPT . EO. 3)  GO  TO  BIO 

AFFECT  A  TIME  SHIFT  IN  SP E 1  CHE R - BRODE  HISTORY 
TO  ALLOW  THE  OVERLAY  TO  BE  PROPERLY  PHASED 
TSHFT  -  (TA-W13/1000  1-TPEB 
DO  800  JT« 1 .NBPTS 

TTIM(JT)  •  TT I  Ml JT  ) - T  SHFT 
800  CONTINUE 

GO  TO  130 

810  CALL  FMAX ( TT I M . NBPT ? . XPMN , XPMX ) 

CALL  FM/Xl PRE SS .NBPTS . YPMN , YPMX ) 

WRITE ( 48 . 640)  XPMN. XPMX . YPMN, YPMX 
840  FORMAT (4E 15. 8) 

IF(IFILT.LT.O)  CO  TO  130 

CALL  for  filters  to  be  executed 

CALL  FL00P ( TTIM. PRESS. DTBN. NBPTS. PF I LT) 

RETURN 

130  WRITE l 48 . 2 10 1  l TTIMl I J  )  .  I J* 1  . NPPTS ) 

WRITEi4e.21C!  ( PRESS ( UI  I  . Ul • 1 .NPpTS  > 

IF  I  IOPT  EC  3  )  GO  TO  350 

IMPULSE 

135  I  T L ( 5  I  *  1 0H  IMPULSE  H 
I T L ( 6  )  =  10HI STORY 
WRITE148.21S)  ITL(5).ITL(6) 

215  F0RMAT12A1O) 

CALL  FMAXITIMP.N1.XIMN.XIMX) 

CAUL  FMAMPIMP.Nl.VIMN.YlMx) 

WRITE (48. 200)  NI.XIMN.XIMX.YIMN.YIMX 
WRITE148.21C)  <TIMP( IV) ,IY«1 .NI ) 

WR I TE( 48.210)  ( P I MP (IT),IT*1,NI ) 

850  I  IMP  -  2 

CALL  I MPULSE ( I  I MP . DTBN , NPPTS . NI ) 

WR I T  E ( 48 , 4  50 )  NI 
IF  ( I0PT.NE.3)  GO  TO  150 
T T L ( 3  )  -  IOHBRODE  IMPU 
I T  L ( 4 )  *  IOHLSE  HISTOH 
I T  L ( 5  )  *  1 0H  Y 

WRITE! 48, 225  )  ITL(3).ITL(4).ITL<5) 

225  FORMAT ( 3A 10 ) 

CALL  FMAXfTIMP.NI .XIMN.XIMX) 

CALL  FMAX ( PIMP .NI , YIMN, YIMX ) 

WRI TE ( 48 . 840 >  X  I MN , X 1  MX , Y IMN . Y I  MX 
150  WRITE ( 48 . 2 10 )  ( T IMP ( KJ I ,KJ* 1 . NI) 

WRI TE ( 48 . 2 10 )  I P I MP (VLI.KL'I.NI) 

IF  (  IOPT  ,NE  .  1  )  GO  TO  175 

FIND  THE  FOURIER  TRANSFORM  ANO  CALCULATE  AMPLITUDE. 

I TL I  5  I  *  1 0H  FOURIER  A 
I T  l ( 6  )  *  10HMPLITUDE  S 
ITlI7)  •  10HPECTPUM 

WPITE(48.225)  ITL(5l.ITL(6).ITL<7) 

CALL  FMAX! FRO. NEF  , XFMN, XFMX  ) 

CALL  FMAX! AMP , NEF . YFMN . YFMX ) 

WR I TE ( 48 . 200 )  NEF.XFMN.XF  MX. YFMN, YFMX 
WRI TE ( 48 . 2 10)  ( FRO!  10  )  , 10*  1  .NEF  ) 

WRITE! 48 . 2 10)  ( AMP! IP ) , IP* 1 ,NEF ) 

175  TOTT  •  DTBN* NBPT S 
C  FREQUENCY  INCREMENT 

DFB  ■  1 . /TOTT 
FOB  ■  0 
WKB  •  0 

NBF  •  NBPTS/2* 1 


1 


DO  349  IK* 1 . NBF 

FRO( LK I  *  0. 

AMP ( LK I  ■  C 

>  r ft ( lk  )  *  o . 

349  CONTINUE 

CALL  FFTRC (PRESS. NBPTS.XFFT, I WKB . WKB ) 

C  AMPLITUDE  SPECTRUM 

DO  50C  KK*  1 , NBF 

FOE  *  F  OE^DFE 
FRO( KK l  *  FOB 

XRB  *  REAL ( XFFT( KK  !  I/NBPTS 
X  I E  •  AIMAGIXF^tIkK)  l/NBPTS 
AMP ( KK  I  *  SORT ( XRB' XRB+X IB*XIB)*TOTT 
500  CONTINUE 
C 

WRITE ( 48 . 450)  NBF 
IF (  I  OPT . NE . 3 )  GO  TO  165 
I  T L ( 3  I  *  10HBRODE  FOUR 
I  T  L  (  4  )  *  10HIER  AMPLIT 
ITL(S)  *  10HUDE  SPECTR 
ITU  6)  =  10HUM 

WRITE(48,235)  1TL(3).ITI(4).ITL(5).ITL(6) 
235  FORMAT ( 4A 10 ) 

CALL  FMAX! FRO . NBF , XFMN, XFMX ) 

CALL  FMAX( AMP . NBF . YFMN. YFMX ) 

WRITE! 4  8 , 840)  XFMN . X  FMX , YFMN. VFMX 
165  WRITE(48.210)  (FRQ< IU) . IU* 1 .NBF > 

WRITE! 48. 210)  (AMP( IE ) . IE* 1 ,NBF ) 

RETURN 

END 

SUBROUTINE  F LOOP ( TT I M. PRESS .DT ,NP . PF I LT ) 


c  this  subroutine  performs  the  looping  required 
C  TO  filter  the  data  or  the  brdde  up  to  seven 

C  TIMES.  FOR  LESS  THAN  SEVEN  FILTER  LEVELS. 

C  FLO  MUST  BE  SET  TO  0.  IN  THE  INPUT  DECK  IN 

C  ORDER  TO  ESCAPE  THE  LOOP. 

C  * . * . 

c 

c 

COMMON  /FJLT/  IFILT,FL0(7). PFQMX ( 7 ) . PFBMX ( 7 ) 
DIMENSION  TTIM!  1  ). PRESS!  1  )  ,PF I LT (  1  ) 

C 

DO  750  JF= 1 .7 

IF ( FLO! JF ) . EQ.O.  )  GO  TO  555 
JFLAG  *  1 

WR I TE ( 48 , 95 )  I F  LAG 

95  FORMAT (15) 

WRITE148.96)  FLO(UF) 

96  FORMAT! F 10. 0) 

DO  725  KF  *  1 , NP 

PFUT(KF  )  ■  0. 

725  CONTINUE 

c 

C  CALL  TO  FILTER 

CALL  FILTER(DT.NP) 

CALL  FMAX(PFILT ,NP . YFMN, YFMX ) 

WRITE! 48, lOO)  YFMN, YFMX 
100  F ORMAT ( 2E 1 5 . 8  ) 

WRITE(4B.  105  I  ( TT IM( LF  I  , LF  =  1  . N& ) 

WRITE ( 48 .  105  )  ( PF ILT (MF  ) . MF  *  1  . NP  ) 

105  FORMAT ( 10E 15.8) 

750  CONTINUE 
555  IFLAG  *  -1 

WR I T  E ( 4  8 , 95 )  IFLAG 

RETURN 

END 


SUBROUTINE  SPLINE (TLAST. NP .TTIM, PRESS ) 


THIS  SUBROUTINE  SETS  UP  A  COSINE  SQUARED  SPLINE 
FUNCTION  AND  APPLIES  IT  TO  THE  FINAL  15%  OF  THE 
PRESSURE  HISTORY  TO  AVOID  A  FREQUENCY  IMPULSE 
IN  TRUNCATED  RECORDS. 


DIMENSION  TT I M( 1  )  , PRE SS (  1 ) 

PIE  =  3.  1415927 
K  =  I  F  I  X (  . 85  *NP ) 

N  *  NP-K-M 
T 1  «  TTIM(K) 

DO  10  «J*  1  .  N 

TFACT  •  (TTIMIK )-T 1  l/< TLAST -T 1  ) 
SFACT  *  COSITFACT-PIE- .5) 

SFACT  ■  SFACT. SFACT 
PRESS(K)  *  PRESS! K  )  «  SF  ACT 
K  *  K*  1 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  I MPULSE ( 1 1 MP , DT . NP , N1 ) 


THIS  SUBROUTINE  CALCULATES  THE  IMPULSE  OF  THE  INPUT 
PRESSURE  DATA  (IIMP  *  1 )  OR  OF  THE  CALCULATED  SPEICHER- 
BRODE  (IIMP  *  2)  EY  SIMPSON'S  APPROXIMATION. 


COMMON  /THIST/  TTIM(6000) ,PRESS(6000) ,TIMP(2999) .PIMP!  2999 ) 
•  PF 1 LT ( 6000 ) 

NTMP  *  NP -  3 
NI  *  NTMP / 2 
DO  90  1*1 .NI 

TIMPI I )  -  O. 

P I MP ( I )  •  0. 

•0  CONTINUE 
IJ  *  O 
SUMIMP  •  O. 

DO  BO  J-3.NTMP.2 

1J  -  IJ+1 

T I MP ( I J  1  •  TT I M( J  ) 

AREA  *  (  PRESS! U- 1  1*4  .  .PRESS(d )*PRESS( 0*1  I >*DT/3. 

SUMIMP  •  SUMIMP+AREA 
D  I  MP  (  I  J  )  ■  SUMIMP 

ac  continue 
return 
END 

SUBROUTINE  FILTER(DT.NP) 

this  SUBROUTINE  filters  the  input  PRESSURE  HISTORY 
(DATA  OR  SPEICHER-ERODE ) .  IT  USES  THE  DICCERENCE 
EOUATIONS  DERIVED  FOR  A  SECOND  ORDER  BUTTERwCRTH 
FIlTER  AS  PRESENTED  Ey  STEARNS.  19’5. 


COMMON  /THIST/  TTIMI 6000) , PRESS(60O0) . TIMP( 2999 ) .PIMP! 2999 ) 
-  PF I LT ( 6000 ) 

COMMON  /COUNT/  I  COUNT .  I  OPT , L F I LT 

COMMON  /FILT  /  IFILT . FLO! 7  )  , PFDMX! 7 ) ,PFBMX( 7 ) 

DATA  LFILT/O/ 

PI  *  3.14  15927 
S2  *  SORT! 2. ) 

LFILT  ■  lfilt-m 


mmmm 
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APPENDIX  D 


LISTING  OF  PROGRAM  FOURPLT 


PROGRAM  FOURPL T ( INPUT .OUTPUT , T APE 5* INPUT . TAPE6*0UTPUT , 

tape9 . plot > 


PROGRAM  FOURPLT  WAS  WRITTEN  TO  PLOT  THE  RESULTS  OF 
PROGRAM  FOURFIT.  IT  READS  INPUT  ON  ASSIGNED  FILE  T  AF  E  22 
AND  PLOTS  SPE1CHER-ERODE  OR  DATA  PRESSURE  AND  IMPULSE 
HISTORIES  AND  FOURIER  AMPLITUDE  SPECTRA  OR  PLOTS 
OVERLAYS  IN  THESE  SAME  DOMAINS  OF  A  DATA  TRACE  AND  ITS 
BEST  FIT  SPEICHER-6R0DE  AS  DETERMINED  BY  FOURFIT. 
WRITTEN  U.  C.  PARTCH  AND  0.  W.  STEEDMAN.  APPLIED 
RESEARCH  ASSOC.,  INC.,  ALBUQUERQUE.  NM ,  FEB  1984. 


COMMON  /ESTIM/  PSOI  ,WI  .PP.W13.PS0F ,WF,F SO 
COMMON  /'PEAK  /  DP.TA.PSO.ALPF 
COMMON  /THI ST /  RSKFT . YS . S , XM 

COMMON  /FILT  /  IFILT.  ITYPE . FLO< 7 ) . FHI ( 7 ) .PFDMXI 7 ) ,PFBMX< 7  ) 

COMMON  /COUNT/  ICOUNT ,  IOPT  . LF  1  LT 
COMMLN  /PLOT V/  1 T L ( 8 ) . 1STLI 8 )  .  IDE 
DIMENSION  X  AR  Y  ( 6000 ) , Y  AR  Y ( 6000  ) 

CALL  GPLOT (  1HU . 7HARAARDS . 7  ) 

CALL  BGNPL ( -  1  ) 

READ! 9 , 100)  IOPT .IFILT 
100  F0RMATI2I51 

IFF IOPT . E0.3)  GO  TO  200 
READ! 9. 120)  <  1 1 L I  I  (.1*1.8) 

120  FORMAT ( BA 1C ) 

IF<  I  OFT  .NE . 2 )  GO  TO  200 

OPTION  2  (DATA  ONLY) 

READ  PRESSURE -T IME  PAIRS 
OR  FILTERED  PRESSURE -TIME  PAIRS 
RE  AO  (  9  .  130)  NEPTS.XPMN.X.PMX.  YPMN.  YPMX 
130  FORMAT* 15 ,4E 15 . 8  ) 

DO  135  IF* 1 .7 

IF(  in  LT  .  LT  .  0  )  GO  TO  ::b 
RE  AD ( 9 .  12  5)  IFlAG 
125  FORMAT!  IS  I 

IFdFLAG.LT. 01  GG  TO  1*2 
PE  AM  9  ’27)  c  L  0  f  IF  ) 

READ!  9  .  172)  YFMfj ,  YPMX 
127  FORMAT ( f 1C  0 ! 

12e  RE  AD  l  9 .  14  0)  (XAPy ( K  I  ,K* 1  .NE^TS ) 

RE  AD 19.  140)  ( Y AP Y ( L  ).L*1  .NEPTS) 

140  FORMAT!  IOC  15.8) 

CALL  PLOTTER l AARY . YARy .NEPTS . XPMN . XPMX , YPMN. YPMX , 1 . 1 . 2 . 1 , 2 ) 
CALL  ENDPL ( -  1 ) 

IF(  IFILT . LT .0)  GO  TO  144 
135  CONTINUE 
112  CALL  GDONE 
STOP  1  1 

READ  IMPULSE-TIME  PAIRS  FOR  OPTION  2 

144  READ(9, 145)  I TL ( 5 ) . I T L < 6 ) 

RE  AD ( 9 . 1 30 )  NE I , X I MN . X I  MX , Y I MN , Y I  MX 
READ (9.1 40 )  ( XARY(M)  ,M> 1 , NE I  ) 

RE  AD (9,1 40  )  ( YARY(N)  . N» 1 . NE I  ) 

145  FORMAT ( 2A 10 ) 

CALL  PLOTTER*  X ARY . Y AR Y . NE I  . X I MN . X I  MX . Y I MN . Y I  MX .  1  .  1  , 3 ,  1 . 3 ) 

CALL  ENDPL  (  -  1  ) 

READ  AMPLITUDE -FREQUENCY  PAIRS  FOR  OPTION  2 
RE  AD ( 9 , 254  )  I T L ( 5 ) . I T L ( 6 ) . I T L ( 7  ) 

RE  AD (9.  130)  NEF.XFMN.XFMX.YFMN.YFMX 
READ! 9 . 140 )  (XARYI  II  ) . 1 1  *  1 . NE  F  ) 

RE AD( 9 , 140 )  (  YARY( JU) . JU* 1 .NEF  ) 

CALL  PLOTTER! XARY . Y ARY . NE F . XFMN . XFMx , YFMN . YFMX . 2 . 5 . 4 , 4 . 3 ) 

CALL  ENDPL ( -  1  ) 

CALL  GOONE 
STOP  777 


71 


OPTION  1  (OVERLAY)  AMO  OPTION  3  ( SPE I  CHER - BRODE ) 

200  READ( 9 , 205 )  PSOF , WF 

RE  ADI  9 , 207 )  DP.TA.RSKFT 
205  FORMAT (2E 15. 8) 

207  FORMAT ( 3E 15 . 8 ) 

I F ( IOPT.NE .  1  )  GO  TO  765 
RE AD( 9 . 1 20 )  ( ISTL(KK)  ,KK=  1 .8  ) 

READ  PRESSURE-TIME  PAIRS  (DATA) 

READ( 9 . 130)  NEPTS.XPMN.XPMX.YPMN.YPMX 
READ (9.1  AO )  ( XARY ( 1T).IT*1.NEPTS) 

READ! 9.  140)  ( YARY ( I W  )  . I W« 1 . NE PT S  ) 

RE  ADI  9 . 766  )  ALPF 
766  FORMAT (  F 10 . 0  ) 

CALL  PLOTTER! XARY , YARY.NEPTS, XPMN . XPMX , YPMN . YPMX . 1 . 1 .2. 1 ,2 ) 

READ  PRESSURE-TIME  PAIRS  OR  FILTERED 
PRESSURE-TIME  PAIRS  ( SPE I  CHE R - BRODE ) 

765  IF< IOPT.EQ. 3)  READI9.120)  I i T L ( I R  ) . I R« 1 . 8 ) 

READ! 9 . 160  )  NPPTS 
160  FORMAT (15) 

IF(I0PT.N£.3)  00  TO  76B 

READ (9,1 70 )  XPMN.XPMX. YPMN, YPMX 

170  FORMAT (4E 15. • ) 

XPMX  •  XPMX/ 4 . 

NBPTS  •  NPPTS/4 

768  DO  234  MF-1,7 

I F ( IOPT . EQ ■ 1  .OR.  IFILT.LT. 0)  GO  TO  171 
REA0(9. 125)  JFLAG 
I F ( JFLAG.LT. 0)  GO  TO  236 
READ( 9 , 127)  FLO(MF) 

RE AD ( 9 ,172)  YPMN. YPMX 
172  FORMAT ( 2E 15 . 8 ) 

171  READ( 9 . 140 )  ( XARY ( LL ). LL« 1 . NPPTS ) 

READ ( 9 , 140)  (YARY(MN) .MN-1 .NPPTS) 

PLOT  SPEICHER-BROOE  ONLY 
I F ( IOPT . EO . 3  )  CALL  PLOTTER 

•  (XARY. YARY, NBPTS, XPMN. XPMX, YPMN. YPMX. 1.1.2, 1.2) 
OVERLAY 

I F ( I  OPT . EO . 1 )  CALL  PLOTTER 

•  (XARY , YARY. NPPTS. XPMN. XPMX. YPMN, YPMX, -1 , 1 .2. 1 .2) 
CALL  ENOPL(-I) 

I F ( IOPT . EO . 1  .OR.  IFILT.LT. 0)  GO  TO  264 
234  CONTINUE 
236  CALL  GO ONE 
ST0P22 


IMPULSE 

264  IF( IOPT . EO. 3  )  GO  TO  280 

READ  IMPULSE-TIME  PAIRS  (DATA) 

RE ADI  9 . 145 )  I TL ( 5 ) . I TL ( 6  ) 

RE ADI  9.  130)  NE I . X I MN , X 1  MX , Y 1 MN , Y I MX 
RE AD( 9 , 140 )  (XARY(NN) ,NN= 1 , NE I  ) 

READ ( 9 . 140 )  ( YARY(MN) ,MN* 1 . NE I  ) 

CALL  PLOTTER! XARY. YARY.NE I , X 1 MN . X I MX , Y I MN , Y I MX , 1 . 1 . 3 . 1 . 3 ) 

READ  IMPULSE-TIME  PAIRS  ( SPE I CHlR- BRODE ) 

280  RE ADI  9 , 160)  N1PTS 

IF ( , OPT . EO . 3 )  READ! 9. 254 )  I TL ( 3  ) . ITL ( 4  ) . ITL ( 5  ) 

254  FORMAT ( 3A 1 0 ) 

IFI I0PT.E0.3)  READI9.170)  X I MN , X I  MX , Y I MN , Y I MX 
READ! 9. 140)  ( XARY ( Id) . I d* i .NIPTS  ) 

READ! 9, 140)  ( YARY ( JI  ). JI ■ 1 , NIPTS 1 

PLOT  SPE I  CHER -BRODE  ONLY 
IFI IOPT . EO. 3)  CALL  PLOTTER 

.  (XARY . YARY , NI PTS . X  I MN . X  I  MX . YIMN. YIMX.  1 .  1 . 3.  1 .3) 

OVERLAY 
NIPTS'NIPTS*  1 

IF ( IOPT . EO ,  1  )  CALL  PLOTTER 

.  (XARY. YARY. NIPTS. XIMN.X1MX.YIMN.YIMX.-1. 1.3.1.31 

CALL  ENOPL  < -  1 ) 


noon  o  °  ^  n  noon  on  ooooo 


FOURIER  AMPLITUDE 
1F( I  OPT . EO. 3  )  GO  TO  340 

RE 40  AMPLITUDE-FREQUENCY  PAIRS  (DATA) 

RE  ADI  9 . 254  )  ITLI5). JTL(6).ITL(7) 

READ! 9 . 130)  NE F . XFMN. XFMX , YFMN. YFMX 
RE  AD  I  9 ,  140)  I XARY ( KJ  I  ,KJ  = 1 ,NEF  ) 

RE  ADI  9 .  140  I  <  YARN  <  JK ) , JK= 1 , NEF ) 

CALL  PLOTTER ( XARt , Y AR Y . NE F , XF MN . XFM/ . Y FMN . Y FMX . 2 . 5 . 4 , 4 , 3 ) 

READ  AMPLITUDE-FREOUENCy  PAIRS  ( SPE 1  CHE R - BR0DE  ) 

340  READI 9 . 160 )  NEF 

IFIIOPT  E0.3)  READ (9.256  )  I TL < 3 > . I TL I  4 1 . I TL ( 5 ) . J TL I  6 > 

256  F ORMAT I  4 A  1 0 ) 

IF! 10PT.E0-3)  READ i  9  . 1701  XFMN . XFMX . YF MN . Y FMX 
RE  AD  19.1 40 )  (XARyIKL  I . K  L  = 1 . NB  F  I 
READI 9 . 140)  (VARy(LKI.LK=1, NEF ) 

SPEICHER-ERODE  ONLY 
IFI10FT.E0.3l  CALL  PLOTTER 

-  (XARY. YARY.NBF. XFMN. XFMX. YFMN. YFMX. 2. 5. 4. 4, 3) 

OVERLAY 

IF! IOPT . EO. 1 )  CALL  PLOTTER 

•  (XARY. YARY.NBF. XFMN, XFMX. YFMN, YFMX. -1.5, 4, 4. 3) 

CALL  ENDPLI -  1  ) 

CALL  GDONE 
END 

SUBROUTINE  PLOTTER ( XARY , Y ARY .NP . XMN.XMX , YMN , YMX . KIND , LBLX . LBLY . 

•  NITSX.NITSY) 

COMMON  /ESTIM/  PSOI ,W1 .PP.W13.PS0F.WF .FSO 
COMMON  /PEAK  /  DP , TA . PSO . ALPF 

COMMON  /FILT  /  IF1LT . 1TYPE ,FL0(7  )  ,FHI I  7 ) .PFDMXI 7  )  .PFBMM 7 ) 
COMMON  /THIST/  RSKFT  .YS.S.X.M 
COMMON  /COUNT  '  I  COUNT . I OPT . LF 3 L T 

COMMON  /PLOT V/  ITL(e ) . ISTL(8 ) . 1DB 

DIMENSION  XARY ( NP ) . Y ARY ( NP ) . LABS (6. 2 ) . LEND (4 . 2 ) . L ABX ( 4 ) . 

•  L  AB  Y I  4  ) 

DATA  ( LABS ( 1 . J ) , J* 1 , 2 )  / 10H  . lOH  TIME  / 

DATA  I LAESI 2 . J ) . J* 1 . 2 )  / 1 0H  .  lOHPRE  SSURE  / 

DATA  I  LABS  I  3 . J  )  . J1 1 . 2 )  / 1 0H  .  10H  IMPULSE  / 

DATA  ( LABS! 4 , J ) . 1 . 2 )  /10H  A . lOHMPL I TUDE  / 

data  1 1  AESI  E  .U  I  .  J*  1 .2  I  /1QH  F  ,  1 0HREQ'JENC  Y  / 

DAI  A  I  l  AE  5  ■  t  .  J  J  .  D  -  '  .  2  •  /-.OP  .  K>H 

data  lendi  * .  /ir.Hi  se:  i 

data  lEn:  ' :  i  ioh1  Mr  /  i  / 

C'ATA  LENM3I  '  10HI  \*-L- SCC  I 

DA 'A  LEND' 4  I  ■  1 0H I  HZ  I 

DATA  LEND!  5  i  MOhi  RADIANS  1  / 

DATA  lfilt/o/ 

WRITE ( 6 . 2300  )  NP . XMN.XMX , YMN. YMX .KIND 
300  F0RMAT(5X.-  ENTERED  PLOTTER  •./. 

•  *NP ,XMN, XMX . YMN, YMX. KIND  «  • ,  I  5 , 4 (  1 X . F7 . 4 ) , I  5 ) 

CALL  H£IGHT(0. 1 ) 

IFIKIND.LT. 0)  GO  TO  200 

DO  10  1*1.2 

LABX(I  )  *  LABSUBLX.I  ) 

LABY ( I  )  *  LABS( LBL Y , I ) 

10  CONTINUE 

L ABX ( 3  )  •  LEND! NIT5X I 
LABY ( 3 )  •  LEND ( N I T  S Y ) 

IFfKIND.E0.2)  GO  TO  100 

•••••  IF  KIND . EO . 1  THEN  PLOT  IS  L I NE AR - L INE AR  ••••• 

50  LINET  •  0 
LINES  •  0 


73 


CALL  SCLMXMN.XMX.  XORG.  XSTP.  XEND) 

CALL  SCL 1 ( YMN. YMX , YORG. YSTP . > END ) 

WRITE (6. 2303 )  XORG . XSTP , XEND . YORG , YST P . Y END 
2303  F0RMAT(2X, -LINEAR  PLOT  • . 6 ( 2X . FB . 4  )  I 

CALL  RLINERI  XORG. XSTP , XEND . YORG, YSTP . YEND , LAEX . LAEY ) 
CALL  DRAWC (XARY.YARY . NP , L I NE T , L I NE S ) 

GO  TO  400 

*•••*  IF  KIND. £0.2  THEN  PLOT  IS  LOG-LOG 

100  LINET  *  0 
LINES  *  0 

CALL  SCL2 ( XMN, XMX , XORG, XCYC .KIND ) 

1 F (KIND . EO. 1 )  GO  TO  50 

CALL  SCL2I YMN. YMX, YORG. YCYC, KIND) 

IF(KIND.EO. 1 )  GO  TO  50 
WRITE ( 6 , 2305  )  XORG. XCYC . YORG. YCYC 
2305  FORMAT ( 5X .• LOG-LOG  PLOT  *  . 4 ( 2X . F B . 4  )  ) 

CALL  LOGLLLl XORG, XCYC . YORG. YCYC, LAEX. LABY) 

CALL  DRAWC ( X ARY . Y ARY, NP. LINET, LINES) 

GO  TO  400 

IF  KIND.LT.O  THEN  PLOT  AN  OVERLAY  ..... 

200  LINET  «  LINET+1 
WR 1 T  E ( 6 , 2307 ) 

2307  FORMAT ( 5X . •  OVERLAY  PLOT  •) 

CALL  BLOFF(IDB) 

CALL  MESSAGt ISTL .80.0  0.6  25 ) 

CALL  MESSAG ( 4 HD AT  A. 4,45.58) 

CALL  STRTPT ( 5 . 2 . 5 . 8 ) 

CALL  CONNPT (5.8.58) 

CALL  MESSAGI 4HF I T  ,4. 4. 5. 5. 6) 

CALL  DASH 

CALL  STRTPT ( 5 . 2 . 5 . 6 ) 

CALL  CONNPT ( 5. e. 5. 6) 

CALL  RESET ( 4HDASH ) 

CALL  DRAWC< XARY . YARv , NP , L INET . L INE S ) 

GO  TO  900 

400  I F ( IOPT . EO . 1  .0«.  IFILT.LT.O)  00  TO  300 
lfilt-lfilt+1 

CALL  MESSAGI 15HL0W  PASS  F I  LTER . 15 . 6 . 5 ,  1 . 0 ) 

CALL  MESSAG( 15HFCUT0FF  (H2)  =  .15.6.5,0.75) 

CALL  REALNOI FLO <  LF1LT I .0.8. 2,0.75) 

300  If < IOPT . EO. 2 )  GO  TO  900 

C A _ L  MESSAGI 13HYIEL0  (KT)  =  .13.6.5.5.5) 

CALL  RE ALNOI WF , 2 ,8 . 2 . f . 5  ) 

Ca„L  MESSAG' 13HPS0  (MPA)  =  .12.0.5.5.26) 

CALL  REALNOI PS0C  . 2 . E . 2 . 5 . 25  ) 

CALL  MESSAGI 13HRANGE  (  KM  I  =  .1". €.6.6.0 

RRP  »  RSKFT.O, 3045 • (WF 23223' ) 

CALL  REALNOI RRR.5.8.2.6.0) 

CALL  MESSAGI 19HP0S .  PHASE  (SEC)  -  .19.6.5.4.751 

DPP  *  DP • 0 . 00 1  *  I WF • • 0 . 3332232  ) 

CALL  REALNOI DPP . 5 . B . 2 . 4 . 75  I 

CALL  MESSAGI 13HT0A  (SEC)  *  .13,6.5,4.5) 

T A A  *  T A • 0  00 1 • ( W F • • 0  3333333) 

CALL  REALNOITAA ,5,8 . 2 . 4 . 5  ) 

IF! IOPT .NE . 1 )  GO  TO  900 
WRITE  I  6 . 666  )  IOPT 
666  F ORMAT ( 2X , • ♦♦♦ I  OPT  * • , 1 5  ) 

CALL  MESSAG(20riLOW  PASS  FID  (HZ)  »  .20.6.5,4.25) 

CALL  REAL  NO (ALPF.0,8.2,4.25) 


#00  CONTINUE 
RETURN 


c 


SUBROUTINE  FMAXI ARY .NA . XMN.XMX ) 
DIMENSION  ARY(NA) 


P 


►  V 


L'-v 

P 

N  „  N 


v.\ 


2300 


10 


WR I T  E ( 6 , 2300  > 

FORMAT ( 5X .• SUBROUTINE  FMAX  * ) 

XMN  =  ARY  I  1  ) 

XMX  *  ARY ( 1 ) 

IF (NA . EQ. 1 )  RETURN 
DO  10  1=2. NA 

IF< XMN. GT . ARY!  1  ) )  XMN  =  ARY  I  I  > 
I F ( XMX , LT . ARY<  I) )  XMX  =  A  R  Y (  I  ) 
CONTINUE 


RETURN 

END 

SUBROUTINE  SCL  1  ( XMN. XMX . AORG. ASTP , AMAX ) 
DIMENSION  S  (  7 ) 


FIND  LINEAR  SCALES  •••> 


2300 


WRITE (6. 2300)  XMN. XMX 
FORMAT (5X. -SUBROUTINE  SCL 1 
SMIN  =  0.00006 
0.00012 
00018 
00024 
00030 
00036 
00060 
00120 


XMN. XMX  =  *.2(F8  4.2X1) 


S(  1  ) 
S<  2  ) 
SI  3  ) 
SI  4  ) 
S(  5  ) 
S  (  6  > 
SI  7  ) 


0. 

0. 

o. 

0. 

0. 

o. 


10 

20 


DIF  =  XMX  -  XMN 
IF  <  DI F . LT . S<  1  ) )  GO  TO  90 
CONTINUE 
DO  10  1=1,7 
IU  =  1 

IFIDIF . LT.S(II)  GO  TO  30 
DO  20  J= 1 ,7 

S(J)  =  S< Jl-10.0 
I F ( S(  1  )  GT . 1 .OE 15)  STOP  1 1  1 
GO  TO  5 


30 


DMA  X  =  S(IU) 
DSTP  =  DMAX/6  0 


DETERMINE  OFFSET 


35 


IF ( XMN  LT .0.0)  GO  TO  60 
DORG  =  0.0 

IF< XMN.LT.DS1P>  GO  TO  99 

01‘SET  =  DSTP 

OFFSET  =  OFFSET*DSTP 

IF ( XMN. GT  OFFSET  )  GO  TO  35 

DORG  *  0FFSET-D5TP 

DMA  X  =  DMAX+DORG 

GO  TO  99 


60 

65 


offset 
OFFSET 
I r  < xm\ 
DORS  = 
DMA/  - 
I  F  i  xv 
I c  lu 


=  0.0 

=  OFFSET-DSTP 
.  L  T .OFFSET  )  GO 
OF' SET 
DMAX -DORG 
LT , DMA> |  GC 
l  T  .  7  I  DMA  )  = 


TO  €  ~ 


1 1 I  I L  EC  ‘  I 
DSTP  -  DMA). 
GO  TU  50 


DMA  • 
l  .  0 


TO  99 
SI  ID*  i  ' 

S  f  x  -  1 C  - .  0 


DIFFERENCE  IS  ZERO 


90 


CONTINUE 
DORG  •  XMN- SMIN 
DMAX  •  XMN+SM1N 
DSTP  •  SMIN/ 3 . 0 


75 


v;v. y; 

,  j.  •  »  »  V  •  U  »  .  M  ,  •  •  ' 


v  v'v'v  *>'*.  /v  v'v  *\*  v%.*%  '*  •-  / 

aji  V  mV  »\k.AWVyV  a  r  ^ 


99  AORG  ■  DORG 
ASTP  •  DSTP 
A  MAX  *  DMAX 

MI  TE  ( 6  , 2303 )  DORG. DSTP .DMAX 
2303  FORMAT ( 5X . •  LEAVING  SCL 1  • . 3( FB . 4 . 2X ) ) 

RETURN 

END 

SUBROUTINE  SCL2(XMN,XMX. AORG. ACVC.KIND) 

SCALE  FOR  LOG-LOG  PLOTS 
MR  I T  E ( 6 , 2  300 )  XMN.XMX 

2300  FORMAT (5X, -ENTER  SCL2  • . 2( F8 . 4 . 2X  )  ) 

I F ( XMN . LT . 1 . OE - B )  GO  TO  80 
1F(XMX.LT . 1 .OE-B)  GO  TO  81 

SMN  «  ALOGIO(XMN) 

SMX  «  ALOGIO(XMX) 

MN  -  I  FIX (SMN) 

IF (SMN. LT . 0. 0)  MN=MN- 1 
MX  *  lFIX(SMX) 

AORG  *  10. • »MN 
DIF  «  (MX-MN)-H 

IF(MN.LT.O  .AND.  MX.LE.O)  DIF  *  MX-MN 
ACYC  *  ABS( 6 . O/DI F ) 

GO  TO  90 

80  WR I TE ( 6 . lOOO )  XMN 

1000  FORMAT  ( SX .  •  XMN  ■=  *.E12.5.»  A  LINEAR  PLOT  WILL  BE  MADE.-) 

GO  TO  82 

81  WR I TE ( 6 . lOOl )  XMX 

1001  FORMAT (5X. -XMX  -  »,E12.5.-  A  LINEAR  PLOT  WILL  BE  MADE.*) 

82  KIND  =  1 

90  CONTINUE 

WRITE (6. 2303)  MN.MX.DIF .AORG.ACVC 

2303  FORMAT (5X. -LEAVING  SCL2  MN . MX . DI F . AORG . ACVC « . 2 1 5 . 3 ( 1 X . F8 . 4 ) ) 

RETURN 
END 

SUBROUTINE  DRAWC ( X . V . NP . L INET . LINES ) 

DIMENSION  X ( NP ) , Y ( NP ) 

WR I TE ( 6 . 2300 )  NP . LINET , LINES 
2300  F0RMAT(5X. -ENTER  DRAWC  NP . LINET . LINES  *  -.315) 

JF(LINET.LE.O)  GO  TO  10 
I F ( LINET . EO . 1 )  CALL  DASH 
I F ( LINET . EO . 2 )  CALL  CHNDOT 
IF(LINET.E0.3)  CALL  CHNDSH 
I F ( LINET . EO . 4 )  CALL  DOT 

10  CALL  CURVE(X.Y,NP, LINES) 

IF(LINET.LE.O)  GO  TO  99 
CALL  RESET! 3HALL  ) 

CALL  HE  I GHT ( 0 .  1 ) 

99  CONTINUE 
RETURN 


SUBROUTINE  RLIN£R( XORG. XSTP . XENO . TORG . YSTP . TEND . LABX . LABY ) 
COMMON  /COUNT/  I  COUNT , IOPT , LF 1 LT 
COMMON  /PLOT V/  ITL(S) . ISTL(B).IDB 
DIMENSION  LABX( 3 ) , LABY ( 3  ) 

C 

WR  j  T  E  <  6 , 2  300  I 

2300  FORMAT  I  !>X  .*  ENTERED  RLINFR . -  ) 

CALL  RAGE  (  1C'.  I .  t  .  5  ) 

CALL  PHYSOPt i.O. 1.01 
CALL  XNAME I LAE>  . 30  I 
CALL  VNAME I LAEY . 30  I 
CALL  AREA2016. 0.6.0) 

IF  (  IOPT  .  EO.  1  I  CALL  BLREC(4.4.5.5. 16,0.6. 1.0) 

IF( IOPT . EO. 1 )  CALL  BLKEY(IDB) 

CALL  MESS AG ( ITL.80.0-0.6.5) 

CALL  GRAF (XORG. XSTP. XEND, YORG. YSTP. YEND) 

CALL  DOT 
CALL  GRID( 1.1) 

CALL  RE  SET ( 3HD0T ) 

C 

RETURN 

END 

SUBROUTINE  LOGLLL ( XOR . XCY . YOR . YCY . LABX . LABY ) 

COMMON  /COUNT/  ICOUNT , I  OPT . LF I LT 
COMMON  /PLOTV/  ITL(8  )  . ISTL(B  ) . IDE 
DIMENSION  LABX ( 3 ) . LABY ( 3 ) 

C 

WRITE ( 6 . 2300 ) 

2300  FORMAT (5X,  'ENTERED  LOGLLL . •) 

CALL  PAGE ( 10.5.8.5) 

CALL  PHYS0R( 1.0.10) 

CALL  XNAME (LABX. 30) 

CALL  YNAME ( LABY . 30 ) 

CALL  AREA2D(6. 0.6.0) 

IF ( IOPT . EO.  1  )  CALL  BLRECI4 . 4 .5. 5. 1 .6.0. 5. 1 .0) 

IF ( IOPT . EO. 1 )  CALL  BLKEY(IDB) 

CALL  MESSAGf ITL.80.0.0.6.5) 

CYC  *  XCY 

IF(YCY.LT.XCY)  CYC  «  YCY 
CALL  LOGLOG( XOR. CYC. YOR. CYC ) 

CALL  DOT 
CALL  GR 1 D ( 1.1) 

CALL  RESET ( 3HD0T ) 

C 

RETURN 

END 


Pso(MPa) 

20 

10 


Amplitude  (MPa -sec) 


- Data 

- Brode 


0.04  008  0.12  0.16 

Time  (sec) 


Pressure  history  for  HEST  record  compared  with  final  fit 
Pso  =  2.95  MPa,  W  =  5.05  KT. 


MPa) 


NORMALIZED  SPEICHER-BRODE  PRESSURE  HISTORIES 


Speicher-Brode  pressure  histories 


NORMflL I  ZED  SPEICHER-BROOE  IMPULSE  HISTORIES 


(SEC/KT 


NORMALIZED  SPEICHER-BRODE  FOURIER  RMPLITUDE  SPECTRR 
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l 
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PEAK  OVERPRESSURE .MPA  = 
NUCLEAR  YIELD. KT  = 

RANGE  FROM  GZ.KM  = 

TIME  OF  ARRIVAL. SEC  = 

POSITIVE  PHASE  DURATION. SEC  = 


. 39600E+02 
. 87 149E+00 
. 24767E-01 
.  15382E-02 
.  14308E  +  00 


+++  LOW  PASS  FIDELITY  (HZ)  «  lOOO.  +  +  + 


figure  16.  Example  FOURFIT  output  for  IOPT  =  1: 

automated  fit  to  0.35  KBAR  DISC  HEST  record 
AB-5  {Speicher-Brode  parameters  listed  on 
file  OUTPUT). 
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HIGH  PRSS  FILTER 
FCUTOFF  I HZ)  -  1QOO 


DISTRIBUTION  LIST 


Defense  Intelligence  Agency 
ATTN:  DB-4C 
ATTN:  D8-4N 
ATTN:  DT-1C 
ATTN:  DT-2 

ATTN:  RTS-2A,  Tech  Library 
ATTN:  RTS-2B 

Defense  Nuclear  Agency 
2  cy  ATTN:  SPSS 
4  cy  ATTN:  STTI-CA 

Defense  Technical  Information  Center 
12  cy  ATTN:  DD 

Dep  of  Defense  Explo  Safety  Board 
ATTN:  Chairman 

Field  Command,  DNA,  Det  2 
Lawrence  Livermore  National  Lab 
ATTN:  FC-1 


Field  Command,  Defense  Nuclear  Agency 


ATTN 

FCPR 

ATTN 

FCT 

ATTN 

FCTT,  W.  Summa 

ATTN 

FCTXE 

Joint  Strat  Tgt  Planning  Staff 

ATTN 

JLK,  DNA  Rep 

ATTN 

JLKS 

ATTN 

JPPFN 

ATTN 

JPST 

ATTN 

JPTM 

Under  Secy  Of  Def  for  Rsch  &  Engrg 

ATTN 

St-at  &  Space  Sys  (OS) 

ATTN 

:  urat  &  Theater  Nuc  For,  B.  Stephan 

DEPARTMENT 

OF  THE  ARMY 

BMD  Advanced  Technology  Center 

ATTN 

:  ATC-T 

ATTN 

•  1CRDABH-X 

BMD  Sys  tens 

Command 

ATTN 

BMDSC-LEE,  R.  Bradshaw 

ATTN 

:  BMDSC-LEE,  R.  Webb 

Dep  Ch  of  Staff  for  Ops  &  Pins 

ATTN 

:  DAMO-NC,  Nuc  Chem  Dir 

Harry  Diamond  Laboratories 

ATTN 

DELHD-NW-P,  20240 

ATTN 

:  DELHD-TA-L ,  81100,  Tech  Library 

US  Army  Armament  Material  Readiness  Conmand 

ATTN 

:  MA  Library 

US  Army  Ballistic  Research  Lab 

ATTN 

:  DRDAR-BLT ,  J.  Keefer 

2  cy  ATTN 

:  DRDAR-BLA-S,  Tech  Library 

US  Army  Concepts  Analysis  Agency 

ATTN:  CSSA-ADL ,  Tech  Library 

US  Ari,.y  Corps  of  Engineers 
ATTN:  DAEN-ECE-T 
ATTN:  DAEN-RDL 

US  Army  Engineer  Div  Huntsville 
ATTN:  HNDED-SR 
ATTN:  HNDED-FO 

US  Army  Engineer  Div  Ohio  River 

ATTN:  ORDAS-L,  Tech  Library 

US  Army  Engr  Waterways  Exper  Station 
ATTN:  J.  Strange 
ATTN:  J.  Zelasko 
ATTN:  Library 
ATTN:  WESSD,  J.  Jackson 
ATTN:  WESSE 

US  Army  Foreign  Science  &  Tech  Ctr 
ATTN:  DRXST-SD 

US  Army  Material  &  Mechanics  Rsch  Ctr 
ATTN:  DRXMR,  J.  Mescall 
ATTN:  Tech  Library 

US  Army  Material  Comnand 

ATTN:  DRXAM-TL ,  Tech  Library 

US  Army  Mobility  Equip  R&D  Cmd 

ATTN:  DRDME-WC ,  Tech  Library 

US  Army  Nuclear  &  Chemical  Agency 
ATTN:  Library 

US  Army  War  College 
ATTN:  Library 

US  Army  White  Sands  Missile  Range 

ATTN:  STEWS-TE-N,  K.  Cunnings 

USA  Missile  Command 

ATTN:  Documents  Section 
ATTN:  DRSMI-RH 

DEPARTMENT  OF  THE  NAVY 

David  Taylor  Naval  Ship  R&D  Ctr 

ATTN:  Tech  Info  Ctr,  Code  522.1 

Naval  Civil  Engineering  Laboratory 
ATTN:  Code  L51,  J.  Crawford 

Naval  Electronic  Systems  Command 
ATTN:  PME  117-21 

Naval  Ocean  Systems  Center 

ATTN:  Code  013,  E.  Cooper 
ATTN:  Code  4471,  Tech  Library 


DEPARTMENT  OF  THE  NAVY  (Continued) 

Naval  Facilities  Engineering  Comnand 
ATTN:  Code  04B 

Naval  Material  Command 
ATTN:  MAT  Q8T-22 

Naval  Postgraduate  School 

ATTN:  Code  1^24  Library 

Naval  Research  Laboratory 

ATTN:  Code  2627,  Tech  Library 
ATTN:  Code  4040,  D.  Book 
ATTN:  Code  4040,  0.  Boris 

Naval  Sea  Systems  Comnand 
ATTN:  SEA-0351 
ATTN:  SEA-09G53,  Library 

Naval  Surface  Weapons  Center 
ATTN:  Code  F31 

Naval  Surface  Weapons  Center 

ATTN:  Tech  Library  &  Info  Svcs  Br 

Naval  War  Col  lege 

ATTN:  Code  E-ll,  Tech  Service 

Naval  Weapons  Center 

ATTN:  Code  266,  C.  Austin 

ATTN:  Code  343,  FKA6A2,  Tech  Svcs 

Naval  Weapons  Evaluation  Facility 
ATTN:  Code  10,  Tech  Library 
ATTN:  R.  Hughes 

Ofc  of  the  Dep  Ch  of  Naval  Ops 
ATTN:  NOP  03EG 

ATTN:  NOP  981 

Office  of  Naval  Research 

ATTN:  Code  474,  N.  Perrone 

Strategic  Systems  Project  Office 
ATTN:  NSP-272 

ATTN:  NSP-43,  Tech  Library 

DEPARTMENT  OF  THE  AIR  FORCE 

Air  Force  Systems  Comnand 
ATTN:  DEB 


Ballistic  Missile  Office/DAA 
ATTN:  ENBF,  D.  Gage 
ATTN:  ENS 

ATTN:  ENS,  W.  Weisinger 
ATTN:  ENSN 

ATTN:  MGEN  A.  Schenker 
ATTN:  PP 


Dep  Ch  of  Staff,  Logistics  &  Engineering 
ATTN:  LEEE 


Dep  Ch  of  Staff,  Rsch,  Dev,  &  Acq 
ATTN:  AF/RDQI 

Foreign  Technology  Division 
ATTN :  NIIS  L i b ra ry 

Rome  Air  Development  Center 
ATTN:  TSLD 

Sacramento  Air  Logistics  Center 
ATTN:  MMEAE,  R.  Dal  linger 

Strategic  Air  Command 
ATTN:  DOWE 
ATTN:  INAO 
ATTN:  NRI/STINFO 
ATTN:  XPFS 
ATTN:  XPQ 

DEPARTMENT  OF  ENERGY 

Department  of  Energy 

Albuquerque  Operations  Office 
ATTN:  CT1D 
ATTN:  R.  Jones 

Department  of  Energy 

Office  of  Military  Application,  GTN 
ATTN:  OMA/RD&T 

Department  of  Energy 

Nevada  Operations  Office 

ATTN:  Doc  Con  for  Tech  Library 

OTHER  GOVERNMENT  AGENCIES 

Central  Intelligence  Agency 
ATTN:  OSWR/NED 


Air  Force  Geophysics  Laboratory 
ATTN:  LWH,  H.  Ossing 

Air  Force  Institute  of  Technology 
ATTN:  ENA 
ATTN:  Library 

Air  Force  Systems  Command 
ATTN:  DLW 

Air  Force  Weapons  Laboratory 
ATTN:  OEX 

ATTN:  NTE,  M.  Plamondon 
ATTN:  NTED,  R.  Matalucci 
ATTN:  NTES-C,  C.  Henny 
ATTN:  SUL 

Air  University  Library 
ATTN:  AUL-LSE 

Asst  Ch  of  Staff,  Intelligence 
ATTN:  IN 


Department  of  the  Interior,  Bureau  of  Mines 
ATTN:  Tech  Library 

Federal  Emergency  Management  Agency 

ATTN:  Asst  Assoc  Dir  for  Rsch,  J.  Kerr 
ATTN:  W.  Chipman,  NP-CP 


NATO  School ,  SHAPE 

ATTN:  US  Documents  Officer 


DEPARTMENT  OF  ENERGY  CONTRACTORS 


University  of 

Ca  1  i  fornic 

Lawrence  Live 

nnore 

National  Lab 

ATTN: 

L-203 

,  R. 

Schock 

ATTN: 

L-203 

,  T. 

Butkovich 

ATTN: 

P.  Coyle 

ATTN: 

Tech 

Info 

Dept  Library 

DEPARTMENT  OF  ENERGY  CONTRACTORS  (Continued) 

Los  Alamos  National  Laboratory 

ATTN:  MS  P364  Reports  Library 
ATTN:  MS530,  G.  Spillman 
ATTN:  R.  Whitaker 

Oak  Ridge  National  Laboratory 

ATTN:  Central  Rsch  Library 
ATTN:  Civil  Def  Res  Proj 

Sandia  National  Laboratories 

ATTN:  Library  &  Security  Classification  Div 

Sandia  National  Laboratories 


ATTN: 

Div  7111,  B. 

Vortman 

ATTN: 

J.  Bannister 

ATTN: 

Org  7111,  L. 

Hill 

ATTN: 

Tech  Library 

3141 

ATTN: 

W.  Roherty 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

Aerospace  Corp 

ATTN:  H.  Mirels 

ATTN:  Library  Acquisition  Ml/199 

Agbabian  Associates 

ATTN:  M.  Agbabian 

Analytic  Services,  Inc,  ANSER 
ATTN:  G.  Hesselbacher 

Applied  Research  Associates,  Inc 


ATTN 
2  cy  ATTN 
2  cy  ATTN 


J.  Bratton 
D.  Steedman 
J.  Partch 


Applied  Research  Associates,  Inc 
ATTN:  D.  Piepenburg 

Applied  Research  Associates,  Inc 
ATTN:  R.  Frank 

Applied  Theory,  Inc 

ATTN:  J.  Trulio 

Artec  Associates,  Inc 
ATTN:  S.  Gill 

AVCO  Systems  Division 

ATTN:  Library  A830 


BDM  Corp 

ATTN 

ATTN 

ATTN 


A.  Lavagnino 
Corporate  Library 
T.  Neighbors 


BDM  Corp 

ATTN:  F.  Leech 
ATTN:  R.  Hensley 

Boeing  Co 

ATTN:  Aerospace  Library 
ATTN:  M/S  42/37,  R.  Carlson 

Boeing  Co 

ATTN:  MS-85-20,  D.  Choate 

California  Research  &  Technology,  Inc 
ATTN:  F.  Sauer 


DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 

California  Research  &  Technology,  Inc 
ATTN:  K.  Kreyenhagen 
ATTN:  Library 

California  Research  &  Technology,  Inc 
ATTN:  Tech  Library 

Cal  span  Corp 

ATTN:  Library 

Carpenter  Research  Corp 

ATTN:  H.  Carpenter 

University  of  Denver 

ATTN:  Sec  Officer  for  J.  Wisotski 

EG&G  Wash  Analytical  Svcs  Ctr,  Inc 
ATTN:  Library 

Electro-Mech  Systems,  Inc 
ATTN:  R.  Shunk 

General  Research  Corp 

ATTN:  Tech  Info  Office 

H-Tech  Labs,  Inc 

ATTN:  B.  Hartenbaum 

Horizons  Technology,  Inc 
ATTN:  R.  Kruger 

I  IT  Research  Institute 

ATTN:  Documents  Library 
ATTN:  M.  Johnson 
ATTN:  R.  Welch 

Information  Science,  Inc 
ATTN:  W.  Dudziak 

Institute  for  Defense  Analyses 
ATTN:  Classified  Library 

Kaman  Sciences  Corp 
ATTN:  Library 
ATTN:  N.  Hobbs 
ATTN:  R.  Ruetenik 
ATTN:  S.  Criscione 

Kaman  Sciences  Corp 

ATTN:  F.  Shelton 
ATTN:  Library 

Kaman  Sciences  Corp 

ATTN:  E.  Conrad 

Kaman  Tempo 

ATTN:  DAS  I AC 

Kaman  Tempo 

ATTN:  DAS  1 AC 

Lockheed  Missiles  &  Space  Co,  Inc 

ATTN:  J.  Weisner,  Dept  80-82 

ATTN:  Tech  Info  Ctr  D/COLL,  D/90-11,  B/ 106 

Martin  Marietta  Corp 

ATTN:  G.  Fotieo 


DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


McDonnell  Douglas  Corp 
ATTN:  R.  Halprin 

McDonnell  Douglas  Corp 
ATTN:  M.  Potter 

Merritt  CASES,  Inc 

ATTN:  J.  Merritt 
ATTN:  Library 

University  of  New  Mexico 
ATTN:  N.  Baum 

Nichols  Research  Corp,  Inc 
ATTN:  N.  Byrn 

Pacific-Sierra  Research  Corp 

ATTN:  H.  Brode,  Chairman  SAGE 

Pacifica  Technology 

ATTN:  R.  Allen 
ATTN:  R.  Bjork 

Patel  Enterprise!,  Inc 
ATTN:  M.  Patel 

Physical  Research,  Inc 
ATTN:  W.  Mendes 

Physics  International  Co 
ATTN:  E.  Moore 
ATTN:  L.  Behmiann 


RSO  Associates 


ATTN 

A.  Kuhl 

ATTN 

D.  Simons 

ATTN 

J.  Lewis 

ATTN 

P.  Haas 

ATTN 

Tech  Info  Center 

ATTN 

W.  Wright 

R&D  Associates 

ATTN 

G.  Ganong 

Rand  Corp 

ATTN 

P.  Davis 

Rand  Corp 

ATTN 

B.  Bennett 

S-CUBED 

ATTN 

D.  Grine 

ATTN 

Library 

ATTN 

T.  Riney 

Science  &  Engrg  Associates,  Inc 
ATTN:  B.  Chambers  III 


Science  Appl ications  International  Corp 
ATTN:  G.  Binninger 


Science  Applications  Inti  Corp 
ATTN:  Tech  Library 

Science  Applications  Inti  Corp 


ATTN: 

J.  Cockayne 

ATTN: 

M.  Knasel 

ATTN: 

R.  Sievers 

ATTN: 

W.  Layson 

Science  Applications ,  Inc 
ATTN:  D.  Maxwell 

Southwest  Research  Institute 
ATTN:  A.  Wenzel 
ATTN:  W.  Baker 

SRI  International 

ATTN:  G.  Abrahamson 

Structural  Mechanics  Associates,  Inc 
ATTN:  R.  Kennedy 

Teledyne  Brown  Engineering 
ATTN:  D.  Ormond 
ATTN:  F.  Leopard 
ATTN:  J.  Ravenscraft 

Terra  Tek,  Inc 

ATTN:  A.  Jones 
ATTN:  Library 
ATTN:  S.  Green 

Tetra  Tech,  Inc 

ATTN:  L.  Hwang 

TRW  Electronics  &  Defense  Sector 
ATTN:  Tech  Info  Center 
2  cy  ATTN:  N.  Lipner 

TRW  Electronics  &  Defense  Sector 
ATTN:  E.  Wong 
ATTN:  G.  Hulcher 
ATTN:  P.  Dai 

Universal  Analytics,  Inc 
ATTN:  E.  Field 

Weidlinger  Assoc,  Consulting  Engrg 
ATTN:  T.  Deevy 

Weidlinger  Assoc,  Consulting  Engrg 
ATTN:  I.  Sandler 
ATTN:  M.  Baron 

Weidlinger  Assoc,  Consulting  Engrg 
ATTN:  J.  Isenberg 
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